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ABSTRACT 



Long-duration gamma-ray bursts (GRBs) require an engine capable of driving a 
jet of plasma to ultrarelativistic bulk Lorentz factors of up to several hundred and 
into narrow opening angles of a few degrees. We use global axisymmetric stationary 
solutions of magnetically-dominated (force-free) ultrarelativistic jets to test whether 
the popular magnetic-driving paradigm can generate the required Lorentz factors and 
opening angles. Our global solutions are obtained via time-dependent relativistic ideal 
magnetodynamical numerical simulations which follow the jet from the central engine 
to beyond six orders of magnitude in radius. Our model is primarily motivated by the 
collapsar model, in which a jet is produced by a spinning black hole or neutron star 
and then propagates through a massive stellar envelope. 

We hnd that the size of the presupernova progenitor star and the radial profile of 
pressure inside the star determine the terminal Lorentz factor and opening angle of 
the jet. At the radius where the jet breaks out of the star, our well-motivated fiducial 
model generates a Lorentz factor 7 ~ 400 and a half-opening angle Oj ~ 2°, consistent 
with observations of many long-duration GRBs. Other models with slightly different 
parameters give 7 in the range 100 to 5000 and Oj from 0.1° to 10°, thus reproducing 
the range of properties inferred for GRB jets. A potentially observable feature of some 
of our solutions is that the maximum Poynting flux in the jet is found at ~ Oj with 
the jet power concentrated in a hollow cone, while the maximum in the Lorentz factor 
occurs at an angle substantially smaller than 0j also in a hollow cone. We derive 
approximate analytical formulae for the radial and angular distribution of 7 and the 
radial dependence of 0j. These formulae reproduce the simulation results and allow 
us to predict the outcome of models beyond those simulated. We also briefly discuss 
applications to active galactic nuclei, X-ray binaries, and short-duration GRBs. 

Key words: accretion, accretion discs - black hole physics - galaxies: jets - hydro- 
dynamics - magnetohydrodynamics (MHD) - methods: numerical 



1 INTRODUCTION 

Models of long-duration gamma-ray bursts (GRBs) require 
the ejected plasma to move at ultrarelativi stic speeds i n 
order to avoid the compactness problem (|Piranl l2005t ). 
The Lorentz factor re quired can be as high as F ~ 400 
ijLithwick fc Sarill200ll ). which necessitates a relativistic en- 
gine capable of launching plasma with an enormous amount 
of energy per particle. Achromatic 'jet breaks' in the GRB 
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afterglow imply a finite geometric opening angle 9j ~ a few 
degrees for a typical long-d uration GRB l|Frail et al.l l200ll ; 
|Piranll2005l : IZeh et alj|2006h . Combined with the observed 
fluence and the known distance to the source, this gives a 
typical event energy of ~ 10 51 ergs, comparable to the ki- 
netic energy released in a supernova explosion. 

An ideal engine for producing ultrarelativistic jets 
with small opening angles, low baryon contamina- 
tion, and high total energies is a rapidly rotat- 
ing black hole threaded by a magne tic field and 
accre t ing at a hyper-Edding t on rate ( Naravan et al.l 
Il992l : iLevinson fc Eichlerl 11993 ; iMeszaros fc Reesl I1997T I 
In such a model, the black hole launches an elec- 
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trom agnetically pure jet via the Blandford-Znajek ef- 
fect jBlandford fc Znaiekl Il977l ). More recently, millisec- 
ond magnetars have been seriously considered as an- 
other possible source of magnetically-dominated out- 



flows (|Usovlll992l; iLvutikovl 120061 ; lUzdenskv fc MacFadvenl 
120071 : Bucciantini et all l2007h ~ The standard alternative 
to th is magne t ic-dri vi ng paradigm is neutrino an nihila- 
tion (|Wooslevl 1 19931 ; IMacFadven fc Wooslevl Il999l ). but 
this mechanism probably does not produce sufficient lu- 
minosity to explain most GRBs |Popham et alj Il999l ; 
|Pi Matteo et alJl2002h . 

Rapidly rotating black holes or millisecond magnetars 
are thought to be th e products of core-collapse (|Wooslevl 
Il993l ; IPaczvnskl 19981) or b i nary collisions of compact ob- 
jects (|Naravan et all Il992l . 120011 ). For failed supernovae, 
the black hole or magnetar is surrounded by an accre- 
tion disc whose corona and wind affect the jet structure 
through force-balance between the jet and the surrounding 
gas. In any core-collapse event the jet must penetrate the 
stellar envelop e which can s i gnificantly modify the struc- 



ture of the jet d Wooslevll 19931; IMacFadven fc Woosle 



let 11 v\ 

Alov et all 'hood; INaravan et all l200ll ; IZhang et al 



1999 



2003 



Alov fc Obergaulingerll2007l ). Indeed, as we demonstrate in 



this paper, it is likely the case that the properties of the 
stellar envelope determine the Lorentz factor and opening 
angle of the jet. 

We seek to understand how magnetized rotating com- 
pact objects can launch jets that become sufficiently ultra- 
relativistic and narrow in opening angle to produce long- 
duration GRBs. To achieve this goal we use the relativis- 
tic ideal magnetohydrodynamical (MHD) approximation, 
which is a valid approxim ation for much of the GRB jet 
(e.g., see iMcKinnevi feoO^. The primary difficulty has been 
in obtaining a self-consistent global model of the jet that 
connects the compact object at the center to large distances 
where the observed radiation is produced. In the context of 
the collapsar model, this means we need a model that goes 
all the way from the black hole or neutron star at the cen- 
ter to beyond the outer radius of the Wolf-Rayet progenitor 
star. 

In the past the MHD 
used in numerous analytical 
standing the physics behind 
tion of relativistic jets. The 
tionary force balance 
alytical studies have 



approximation has been 
efforts aimed at under- 
acceleration and collima- 
MHD equations for sta- 
are highly non- linear, and so an- 
been mostly confined to special 



cases for which the equation s can be simplified, e.g . 
for particular field geometries dBlandford fc Znaiekl 1 19771 ; 
IBeskin et ail 1 19981; iBeskin fc Nokhrinal I2006T). for asymp- 



totic solutions (lAppl fc Camenzind 19931; Begelman fc Lil 



1 1994 lLovelace fc Romanoval 120031 ; iFendt fc Quvedl |2004 ). 
or for self-similar solutions that allow variable sepa - 



ration (jContopoulos fc Lovelace _ 19941; Contopoulosl Il995l : 



IVlahakis fc Konigll l2003al lbl; INaravan et all tooll Semi- 
analytical methods using finite element, iterative relaxation, 
or shooting techniques have also been use d to find jet solu- 
tions, such as for s pinning neutron s tars (ICamenzin dl ll987l ; 
lLovelace et~aiX 2006) and black holes (lFendtlll997l ). Such an- 
alytical studies are useful since sometimes one finds families 
of solutions that provide sig nificant insight into the general 
properties of solutions (e.g.. INaravan et al.ll2007l ). 

Time-dependent simulations complement analytical 



studies by allowing one to investigate a few models 
with much less restrictive assumptions. In particular, 
much recent insight into the accretion-jet phenomenon 
has been achieved within the framework of general 
relativistic magnetohydrodynamics (GRMHD) via time- 
dependent numerical simulations (|Pe Villiers et all [20031; 
McKinney fc Gammiel |2004|; |McKinnevll2005bl ; iKomissarovl 



20051 ; I Alov fc Obergaulingerl|20 07i). Indeed, numerical sim- 



ulations of accretion have successfully reproduced colli- 
ma ted relativistic ou tflows with Lorentz factors reaching 
10 (|McKinnevl2006bl ). Within the collapsar model, GRMHD 
simulations show tha t magnetized jets c an be produced 
during core-collapse ( Mizuno et all 2004 ■ Liu et ail 120071 ; 



iBarkov fc Komissarovl 20071 ; Stephens et all 20081 ). However, 
no MHD simulation of core-collapse has yet demonstrated 
the production of an ultrarelativistic jet. Computationally, 
such simulations are prohibitively expensive due to the need 
to resolve vast spatial and temporal scales while at the same 
time modeling all the physics of the black hole, the accretion 
disc, the disc wind, and the stellar envelope. A more prac- 
tical approach, one that we take in this paper, is to replace 
the real problem with a simplified and idealized model and 
to explore this model over the large spatial and temporal 
scales of interest for long-duration GRBs. Such an approach 
will hopefully demonstrate how ultrarelativistic jets can be 
produced and will help us assess the applicability of the 
mechanism to the collapsar model. 

In the present work we obtain global solutions of ultra- 
relativistic magnetically-dominated jets via time-dependent 
numerical MHD simulations in flat space-time (no grav- 
ity). We focus o n the relativistic maqnetodynamic al, or 



1974 



1977 



force-free, regime fGoldreich fc Julian! 19691; Qkamoto] 
Blandfordl Il976l; lLovelace! 1 197d; iBlandford fc Znaiekl I 
MacDonald fc Thornelll982l; IFendt et allll995l ; IKomissarovl 
200ll . 120021 ; iMcKinnev! l2006al ). which corresponds to a 
magnetically-dominated plasma in which particle rest- 
mass and temperature are unimportant and are ig- 
nored. This is a reasonable model for highly magne - 
tized flows l|Blandford fc Znaiekl [19771 ; IMcKinnev! l2006bl) . 
The model parameters we co nsider are motivated 
by p r esupernova stella r models (IMacFadven fc Wooslev 



1999 



2005 



Alov et all l2000l: IZhang et all 120031 ; iHeger eTal 



Zhang et ail [2007) an d GRMHD simulati o ns of 
turbulent accretion d iscs (McKinney fc Gammiel |2004| ; 
IMcKinnev fc Naravanl l2007al lbTK We compare our nu- 
me rical solutions against self-similar solutions derived 
by INaravan et al.l (I2007T ) and obtain simple physically- 
motivated formulae for the variation of the Lorentz factor, 
collimation angle, and Poynting flux along the axis of the 
jet and across the face of the jet. Based upon our simula- 
tions and analytical scalings, we suggest that the terminal 
Lorentz factor of GRB jets is determined by the size and 
radial pressure profile of the progenitor star rather than the 
initial magnetization, for a large range of initial magnetiza- 
tions. 

In [J2] we discuss the problem setup and give a brief 
overview of our numerical method. In 331 we present the 
numerical results and interpret them in terms of analyti- 
cal scalings. In Sj4]we make a comparison to other models. 
In [JS] we discuss astrophysical applications of our models, 
and in f|6]we give a brief conclusion. Readers who are not 
interested in the details may wish to look at Figures [T] - [3] 
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Figure 1. Cartoon of the large-scale structure of a GRB source 
(not to scale). The major elements are a central engine which 
launches a polar magnetically-dominated ultrarelativistic jet, and 
a gaseous stellar envelope (gray shading) which confines the jet. 
The central engine may be an accreting rapidly rotating black hole 
or a millisecond magnetar. For a failed supernova, there could also 
be a disc wind which may additionally confine the jet. 



and to read S}5] In Appendix[A]we introduce an approximate 
model of force-free jets and present a comprehensive discus- 
sion of the analytical properties of these jets. In Appendix [Bl 
we discuss the kinematics of any (dynamically unimportant) 
plasma that may be carried along with a force- free jet. 



2 MOTIVATION, PROBLEM SETUP AND 
NUMERICAL METHOD 

As depicted in Figure [T] a crucial aspect of the collap- 
sar GRB model is that the central engine must produce 
a jet powerful enough to penetrate the stellar envelope. 
The interaction between the stellar envelope and the jet 
is found to be quite complex in time-dependent hydrody- 
namical numerical simulation s of jets injected at an inlet 
within the presupernova core (lAlov et al" 2000; Zhang et al.l 
120031 ; iMorsonv et~al]|2007l ; IWang et al.ll2007l ). These simu- 
lations show that the jet collimates and accelerates as it 
pushes its way through the confining stellar envelope, thus 
suggesting that the envelope plays a crucial role in deter- 
mining the opening angle and Lorentz factor of the flow 
that emerges from the star. If the collapsar system forms an 
accreting black hole, then the ultrarelativistic jet may be ac- 
companied by a moderately relativistic d isc wind that may 
provid e additional collimation for the jet (iMcKinnevl feoOSb . 
l2006bl ). We note that the larger the radius of the progen- 
itor star and/or the denser the stellar envelope, the more 
energy is required for the jet to have to penet rate the stellar 
envelo pe and reach the surface of the star. iBurrows et al.l 
i|2007l ) find that a relativistic jet in the collapsar scenario 
may be preceded by a non-relativistic precursor jet that 



might clear the way for the second, relativistic jet. In the 
magnetar scena rio, the stellar envelope is the primary col- 
limating agent (Uzdcnsky & MacFadycn 2007). Eventually, 
one would like to study individually the collapsar and other 
models of GRBs. However, at the basic level, all models are 
fundamentally similar, since they involve a central magne- 
tized rotating compact object that generates a jet confined 
by some medium (e.g., Fig. [TJ. 

Figure [2] shows our idealized approach to this problem. 
We reduce the various scenarios to a rigidly rotating star of 
unit radius surrounded by a razor-thin differentially rotating 
disc. Magnetic field lines thread both the star and the disc. 
We identify the field lines emerging from the star as the 'jet' 
and the lines from the disc as the 'wind.' Both components 
are treated within the magnetodynamical, or force-free, ap- 
proximation. That is, they are taken to be perfectly con- 
ducting, and we assume that the plasma inertia and thermal 
pressure can be negl ected. In terms of the standard magnet i- 
zation parameter a (|Michelll 19691 : iGoldreich fc Julianlll970h , 
we assume a — *■ oo. In this idealized model, the force-free 
disc wind plays the role of the stellar envelope (plus any 
gaseous disc wind) that collimates the jet in a real GRB 
(Fig. ID. 

In the context of the collapsar picture, the 'wind' region 
of our idealized model can be considered as a freely moving 
pressure boundary for the jet. The jet boundary in our sim- 
ulations is able to self-adjust in response to pressure changes 
within the jet, and thus the boundary is able to act like the 
stellar envelope in a real collapsar. Notice that replacing the 
stellar envelope with our idealized magnetized 'wind' is a 
good approximation because in ideal MHD the wind region 
could be cut out (along the field line that separates the jet 
and wind) and replaced with an isotropic thermal gas pres- 
sure. The problem would be mathematically identical if the 
material in the wind region were slowly moving, as is true 
for the stellar envelope. Even when the material in the wind 
region is rapidly rotating, we show later that the pressure 
in the wind region changes very little, so the approximation 
still remains valid. Since the only importance of the wind 
in our model is to provide pressure support for the jet, we 
adjust our disc wind to match the expected properties of the 
confining medium in a collapsar. 

We work with spherical coordinates (r,9,<p), but we 
also frequently use cylindrical coordinates R = rsin(9, z — 
r cos 8. We work in the units c = ro = 1, where c is the speed 
of light and ro is the radius of the compact object. There- 
fore, the surface of the compact object is always located at 
r = 1, and the unit of time is ro/c. 



2.1 Jet Confinement 

Most of the jet power output from a BH accretion sys- 
tem moves along field lines that originate from the com - 
pact object (|McKinnev fc Gammiell2004lMcKinnev||2006bh . 
i.e., along the zone that we identify as the 'jet' in our 
model (Fig. [5J ■ A crucial factor which determines the degree 
of acceleration and collimation of the jet is the total pressure 
support provided by the surrounding medium through gas 
pressure, magnetic pressure, ram pressure, or other forces. 
We parameterize the initial total confining pressure (at time 
t — 0, which is the starting time of the simulation) as a 
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Figure 2. Idealized model studied in this paper. The thick solid 
lines in the upper panel show an azimuthal cut through a compact 
star surrounded by a razor-thin disc. The star and the disc are 
threaded by magnetic field lines, which are shown as thin solid 
lines. The magnetized plasma above the star and the disc is as- 
sumed to be perfectly conducting and to have an ultra-high mag- 
netization parameter. Arrows show the direction of the poloidal 
electric current. The thick dashed line indicates the field line that 
separates the jet from the disc wind. The disc wind provides pres- 
sure support for the jet and plays the role of the gaseous stellar 
envelope in Fig. [T] The degree of pressure support is adjusted by 
varying the magnetic field strength profile in the disc. The lower 
panel shows the angular frequency of rotation of field lines as a 
function of the cylindrical radius of their foot-points. 



power- law: 



Pamb oc r 



(1) 



Near the compact object the accretion disc it- 
self can provide the jet w ith support. For example, 
IMcKinnev fc Naravanl (|2007al lbh found via GRMHD simula- 
tions of magnetized accreting tori that the wind from their 
torus has an ambient pressure support that approximately 
follows a simple power-law with a ~ 2.5 for two decades 
in radius. At a larger distance from the compact object the 
disc wind is expected to become less effective, and the am- 
bient pressure in the case of a GRB is presumably due to 
thermal and ram pressure of the stellar envelope. A ccord- 
ing t o a simple free-fall model of a collapsing star (|Bethel 
1990), for which density and velocity scale as p oc r~ 3 ^ 2 and 
v oc r _1//2 , the ram pressure varies with radius as ~ r -2 ' 5 , 
identical to the GRMHD disc wind result. Moreover, hy- 
drodynamic simulations of GRB jets show that the internal 
thermal pressure als o has the same sc aling, a ~ 2.5 (see, 
e.g., Model JA-JC in lZhang et aDl2003h . 

In our model, the vertical component of the magnetic 
field at the surface of the disc is taken to vary as a power-law 
with radius, 



This is our boundary condition on the field at z = 0, R > 1. 
If v — 1, the wind has a paraboloidal shape and the magnetic 
pressure has a power-law scaling r~ 2 , whereas if v — 0, the 
wind cor responds to a split monopo l e with pressure varying 
as r~ 4 (iBlandford fc Znaiekl Il977l ; IMcKinnev fc Naravanl 
l2007al lbl; iNaravan et al.ll2007l ). For a general value of v, the 
magnetic pressure in the wind is very close to a power-law, 
r" Q , with 



Since we wish to have a 
model we choose 



a = 2(2- u). (3) 
2.5, therefore for our fiducial 



2 - a/2 = 0.75. 



(4) 



We have considered many other values of u, but focus on 
two other cases: v — 0.6, 1. 



2.2 Model of the Central Compact Object 

We treat the central compact object as a conductor with a 
uniform radial field on its surface, i.e., as a split monopole. 
The compact object and the field lines rotate at a fixed an- 
gular frequency Qo, and it is this magnetized rotation that 
launches and powers the jet. We neglect a l l gravi tational 
effects. As shown in IMcKinnev fc Naravanl (|2007al lbr), this 
is a good approximation since (relativistic) gravitational ef- 
fects do not qualitatively change the field geometry or so- 
lution of the magnetically-dominated jet even close to the 
BH. Also, jet acceleration is known to occur mostly at large 
distances from t he compact object for ele ctromagnetically- 
drivenjets (e.g., iBeskin fc Nokhrinall2006r i. 

For a spinning BH the angular frequency of field lines 
in the magnetosphere is determined by general relativistic 
frame dragging in the vicinity of the hole. This causes the 
field lines threading the BH to rotate with a frequency ap- 
proximately equal to half the rotation frequency of the BH 
horizon, 



fi(a)^0.5n H (a) = -^ = ^, 



(5) 



where the dimensionless Kerr parameter a describes the 
BH spin and can take values between —1 and 1. In our 
chosen units the radius of the BH horizon, rn — (1 + 
~~a?)GM/c 2 , is unity (see 32) ■ Equation (El) is for a 



monopo le field threading the horizon l)McKinnev fc Naravanl 
l2007al lbl). For field geometries other than a monopole, 
the field rotation frequency does not remain exactly 
constant on the BH horizo n (jBlandford fc Znaiekl 11971 
IMcKinnev fc Naravanl l2007bh . For example, for a parabol- 
ical field geometry, Q near the poles is smaller than 
in the monopole case by a factor of two. We do not 
consider this effect in the current study. According to 
IMcKinnev fc Naravanl (|2007bl ) it should not change our re- 
sults significantly. 

Various studies of BH accretion systems sug- 
gest that rapi dly spinning BHs (a > 0.9 ) are 
quite common (Gammie et al.l 120041 : IShafee et al. I 120061 : 



iMcClintock et al 



20061 ). Therefore, for our fiducial model, 



we consider a maximally spinning BH with Kerr parameter 
a = 1, so we choose 



B Z (R) oc R 1 " 



const. 



(2) 



Q{a = 1) = 0.25. 



(6) 
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This is the maximum frequency that field lines threading a 
BH can have in a stationary solution. 

Even though we primarily focus on the case of a maxi- 
mally spinning BH, we also apply our model to magnetars. 
A magnetar with a characteristic spin period of 1 ms and 
a radius of 10 km has spin frequency Ons ~ 0.21 in the 
chosen units (unit of length ro = 10 6 cm and unit of time 
ro/c ~ 3.3 x 10 -5 s), and so is comparable to a maximally 
rotating black hole. Thus, fto = 0.25 is a reasonable approx- 
imation for either rapidly rotating black holes or millisecond 
magnetars. 

2.3 Astrophysical Problem Setup: Models A, B, 
and C 

Since we study axisymmetric magnetic field configurations, 
it is convenient to separate poloidal and toroidal field com- 
ponents, 

B — Bp + B v (p. (7) 

It is further convenient to introduce a magnetic field 
stream fu nction P to describe the axisymmetric poloidal 
field B„ dOkamotdll978l ; iThorne et ail 1 19861 ; iBeskinl 1 19971 ; 
iNaravan et al.ll2007t ). 

1 dP 1 dP ~ 

B v = ^-n%^^-^-e. (8) 

r 2 sin 9 do r sin 9 or 
This representation automatically guarantees V ■ B = 0. 
Here f, 6, and (p are unit vectors in our spherical co- 
ordinate system. The stream function gives the magnetic 
flux cfr enclosed by a toroi dal loop passing through a point 
(r, 6<) (|Naravan et al.ll2007l 1. 

$(r,0) = 2irP(r,6). (9) 

We perform the simulations over the region (l,r max ) x 
(0, tt/2). We initialize the simulation with a purely poloidal 
initial magnetic field, 

P = r"(l - cosS). (10) 

This initial field corresponds to a split monopole field con- 
figuration at the compact object (constant \B r \) and has 
a power-law profile for the vertical component of the field 
on the disc. For our fiducial model A, we take v = 0.75, 
as explained in £12.11 This magnetic field configuration has 
a confining pressure varying as r -2,5 and is approximately 
an equilibrium nonrotating jet solution as we show in Ap- 
pendix [XT] 

We consider both the surface of the compact object, 
r — 1, and the surface of the disc, 9 = tv/2 (z = 0), to 
be ideal conductors. The number of quantities we fix at 
t hese boundaries is consistent with the counting argument 
of lBogovalovl (ll997T l . A paraphrasing of this argument is that 
the number of quantities that one should relax at the bound- 
ary of a perfect conductor equals the number of waves en- 
tering the boundary. In our case of a sub-Alfvenic flow there 
are two waves entering the boundary: an incoming Alfven 
wave and an incoming fast wave. Thus, we leave the two 
components of the magnetic field parallel to the conductor 
unconstrained, and we only set the normal component of the 
field. 

We set the values of two magnetic field drift velocity 
components (perpendicular to the magnetic field) through 



the stationarity condition (|Naravan et al.l l2007i ). For the 
compact object we choose a constant angular velocity ro- 
tation profile flo — 0.25 (eq. |5J), and for the disc we choose 
a Keplerian-like rotation profile, 

n Aisc (R, 2 = 0) = n R~ B , P = 3/2. (11) 

Where the compact object meets the disc, the magnitude 
of the disc angular frequency per unit Keplerian rotation 
frequency is Q/Qk ~ 1/2 as consistent for millisecond mag- 
netars of 10 km size and mass I.4M0 and consistent with 
GRMHD simulations for near an a = 1 black hole (see, e.g., 
iMcKinnev fc Naravanl l2007al lbl). Therefore, near the com- 
pact object our model is more accurate than a precisely 
Keplerian model while keeping Q everywhere continuous. 
We use an antisymmetric boundary condition at the polar 
axis 9 = and an outflow boundary condition at the outer 
boundary r = r max = 10 8 . Since our time-dependent solu- 
tions never reach this artificial outer boundary, our results 
do not depend upon the details of the boundary condition 
there. 

We discuss results obtained with this fiducial model A 
in £13.11 We also discuss two other models that differ from 
model A only by the value of v. model B, which has v — 1 
(parabolic field lines), and model C, which has v = 0.6. We 
have run a number of other models that we mention briefly 
in 3331 

2.4 Numerical Method 

We use a Godunov-type scheme to numericall y solve the 
time-d ependent force-free equations of motion (Mc KinnevI 
2006a). Our code has been successfully used to model 
BH and neutron star magne to spheres (McKinncv 2 006al lcl; 
IMcKinnev fc Naravanll2007al lbl; INaravan et alj|2007l ) . 

To ensure accuracy and to properly resolve the jet, we 
use a numerical grid that approx imately follows the mag- 
netic field lines in the jet solution (|Naravan et al.ll2007T ). We 
are thus able to simulate the jet out to large distances with- 
out making significant errors in the solution. For the three 
models, A, B, C, discussed in this paper we used a resolution 
of 2048x256. Since our grid follows the poloidal field lines, 
the above resolution corresponds to an effective resolution of 
about 2048x100, 000 in spherical polar coordinates. A com- 
parison with lower-resolution runs shows that these models 
are well converged. In particular, over 6 orders of magnitude 
in distance from the compact object, the shape of poloidal 
field lines is accurate to within 89/9 < 10% (see Sj3j) . 

In order to speed up the computations, we use a time- 
stepping technique such that only the non-stationary region 
is evolved. This is achieved by defining the active section, 
where the evolution is performed, to be the exterior to a 
sphere of radius r sta t = fstatrf, where t is the time of the 
simulation, c is the speed of light, and £ s tat = 0.01. We set 
the electromotive forces at all boundaries of the active sec- 
tion to zero. If the initial condition is a force-free solution, 
then this procedure is mathematically justified even within 
the fast critical surface since the time-dependent solution 
only contains outgoing waves, so the solution rapidly settles 
to its final state behind the wave. In all cases our initial 
condition is an exact solution or is close enough to the exact 
solution to avoid significant ingoing waves. We also experi- 
mentally verified that by not evolving the solution interior 
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Figure 3. Poloidal magnetic field lines, shown as solid lines, overlaid on the colour-coded absolute magnitude of the enclosed poloidal 
current \I\ = \RB v /2\ (upper panels), and the colour-coded Lorentz factor 7 (lower panels). The results are for the fiducial model A, 
which corresponds to u = 0.75, Oo = 0.25, /3 = 1.5. Left (right) panels show the near (far) region of the jet solution. The thick dashed 
lines indicate the position of the field line that separates the jet from the disc wind (see Fig. [2}, and the thick solid lines show the position 
of the Alfven surface, QR = 1. Red (blue) colour corresponds to high (low) values of the plotted quantities. Note that, for any distance 
from the compact object, the maximum \I\ is nearly coincident with the jet-wind boundary. However, the maximum Lorentz factor is 
found inside the jet, closer to the axis (see Fig. 1121 . 



to the active section, we make an error of less than one part 
in ten thousand. The use of grid sectioning speeds up the 
simulations by a factor of up to a thousand since it allows 
us to use a larger time step. iKomissarov et all (|2007l ) used 
a similar approach. 



3 SIMULATION RESULTS 

We first present simulation results for our fiducial model A 
and analyze them morphologically. Then, we develop an an- 
alytical model of the jet structure and use it to gain insight 



into jet acceleration and collimation. Finally, we consider 
the other two jet models, B and C, and discuss the variation 
of jet properties with 6. 

3.1 Fiducial model A 

Model A consists of a compact object of unit radius rotat- 
ing with an angular frequency fio = 0.25, surrounded by a 
Keplerian-like disc [ft — 3/2). On the surface of the com- 
pact object, the radial component of the magnetic field is 
uniform, B r = 1, while at the disc, the vertical component 
of the field varies with radius in a self-similar fashion with 
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Figure 4. Radial dependence of the Lorentz factor 7 in the 
fiducial model A for two field lines. One field line starts from 
the compact object at an angle 6{ p 73° (indicated with thick 
lines), and the other starts at 6>f p PS 21° (thin lines). Solid lines 
show the numerical solution, dashed lines show the analytical 
approximation l|13|l with C = \/3 (the solid and dashed lines are 
virtually indistinguishable for 9{ p = 21°), and dotted lines show 
the individual scalings given in (1 2 1 I t and ( 1221 ). Note that the field 
line with 8{ p = 73° accelerates quickly as it moves away from the 
compact object but it then switches to a slower second regime 
of acceleration. In contrast, the field line with 8{ p = 21° begins 
accelerating only after it has moved a considerable distance from 
the compact object. However, it then maintains a rapid rate of 
acceleration without switching to the second acceleration regime. 
When the jet reaches the outer edge of the simulation at r ~ 2 X 
10 6 , this field line has a very large bulk Lorentz factor 7 > 1000, 
whereas the field line with 9{ p = 73° has a smaller 7 ~ 500. Thus, 
the jet develops a fast core surrounded by a slower sheath. 



index v = 0.75, i.e., B Z (R, z = 0) oc R v ~ 2 = R~ 125 (eq.\W$. 
Starting with this purely poloidal initial field configuration, 
we have run the force-free simulation for a time equal to 
10 7 ro/c. At the end of the calculation we obtained a time- 
steady solution out to a distance of 2 x 10 6 ro- We describe 
below the properties of this steady solution. 

The panels in Fig. [3] show the poloidal field structure 
of model A in steady state. The poloidal field in the fi- 
nal rotating state is nearly the same as in the initial non- 
rotating stat e, just as was s e en for the self-similar solutions 
discussed in iNaravan et all l|2007h . This is despite the fact 
the final steady solution has a strong axisymmetric toroidal 
field B v (r,9), which is generated by the rotating boundary 
conditions at the star and the disc. 

The toroidal field at any point is related to the total 
enclosed poloidal current l(r,0) at that point by Ampere's 
Law, 

I(r,6) = RB v {r,6)/2 sC 0, 7? = rsin<9. (12) 
The enclosed current is negative because we have a posi- 




Figure 5. Field line shape (upper panel) and comoving magnetic 
pressure (lower panel) in the fiducial model A for a field line that 
starts from the compact object at 9{ p £3 73°. Solid lines show the 
results from the numerical simulation, and dashed lines show the 
analytic scalings for the non-rotating solution. Field lines with 
foot-points in the disc (wind) show better agreement with the 
analytic scalings. For example, the magnetic pressure along a field 
line with rf p > 500 (not shown) is indistinguishable from the 
dashed line in the lower panel. 



tive B z and positive £1, so that B v is negative. The colour- 
coding in the upper panels of Fig. [3] indicates the absolute 
magnitude of the enclosed current as a function of position. 
As expected, we see that I is constant along field lines, which 
corresponds to RB V being constant. More interestingly, we 
see that at any r, the absolute value of the enclosed cur- 
rent starts at zero, increases as we move away from the axis, 
reaches a maximum value, and then decreases back to zero. 
The maximum in the absolute enclosed current corresponds 
to a transition from a negative current density (inward cur- 
rent) to a positive current density (outward current). This 
transition is coincident with the field line that originates at 
rtp = 1, #f p = 7r/2 and that defines the boundary between 
the 'jet' and the 'wind' in our model (see Fig. [2}. 

As a result of rotation, the solution develops a poloidal 
electric field in the lab-frame (or inertial frame). The elec- 
tric field strength at each point is equal to QRB P , where fi is 
equal to the angular frequency at the foot-point of the local 
field line. The electric field gives an outward Poynting flux 
S — E x B/4n which we discuss later. It also gives a drift 
speed v — E/B, and a corresponding Lorentz factor 7. The 
colour-coding in the lower panels of Fig.[3]indicates the vari- 
ation of 7 with position in the steady solution. The Lorentz 
factor reaches up to a maximum ~ 1000 in this particular 
model. As Fig. [3]shows, the acceleration proceeds gradually 
and occurs over many decades in distance from the compact 
object. 

Note that, at a given distance from the compact object, 
the maximum Lorentz factor is not achieved at either the jet- 
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wind boundary or on the axis, but occurs at an intermediate 
radius inside the jet. For instance, at z — 5 x 10 s , 7 is 
maximum at R ~ 3 x 10 3 , whereas the jet-wind boundary 
is located at R ~ 8 x 10 3 . Thus, the jet consists of a slow 
inner spine, fast edge, and a slow o uter sheath which a c tually 
contains most of the power density. iKomissarov et ail (|2007T ) 
apparently observed this 'anomalous' effect in one of their 
solutions. In the next subsection we explain the origin of the 
effect and quantify it. 

Figure U shows the variation of the Lorentz factor with 
distance along two field lines emerging from the compact 
object. The field line that starts closer to the equator, with 
6f p — 73°, undergoes rapid acceleration once it is beyond a 
distance ~ 10. However, at a distance ~ 10 3 it switches to 
a different and slower mode of acceleration, reaching a final 
7 ~ 500 at r = 2 x 10 6 . In contrast, the field line that starts 
closer to the axis at 9f p = 21° does not begin accelerating 
until a distance ~ 100. It then accelerates rapidly almost 
until it reaches the outer radius (there is a hint of a transition 
to the slower acceleration mode near the end), by which 
point it has a larger Lorentz factor ~ 1000 than the other 
field line. This inverted behaviour is what causes the natural 
development of a fast structured spine and slow sheath that 
contains most of the power density. 

The upper panel in Figure [5] shows, for the steady state 
solution, the variation of polar angle 9 as a function of dis- 
tance along the field line that starts at 9fp = 73° . The dashed 
line shows the corresponding quantity for the initial purely 
poloidal field with which the simulation was started. We see 
that the final field shape is mildly perturbed by rotation. 
However, even at a distance of 2 x 10 6 the change in 9 is no 
more than a factor of 2. 

The lower panel in Figure [5] shows the variation of the 
comoving pressure with distance along the same field line. 
The pressure varies as r~ 2 ' 5 (the dashed line), the desired 
dependence, for distances up to r ~ 10 3 or so. Beyond this 
distance, the pressure variation in the jet becomes a little 
shallower. The change occurs in the region where the slower 
mode of acceleration operates (see Fig. [5} . We explain this 
behaviour in the next subsection. We note that this change 
in pressure inside the jet does not affect the confining pres- 
sure profile of the external medium/wind which stays the 
same as in the initial configuration, and varies as r -2 ' 5 . 

The model A simulation described here is well- 
converged: the angular frequency of field line rotation fl 
and the enclosed poloidal current I axe accurately preserved 
along each field line, as they s hould be in a stationary ax- 
isymmetric force-free solut ion (Mestel 1961; lOkamotoll 19781 ; 
iThorne et ailll986l ; lBeskinlll997i ; lNaravan et al.ll2007f ). Even 
though the simulation domain extends over more than six 
decades in radius, these field-line invariants are conserved to 
better than 12%, see Fig. [6] 



3.2 Comparison to Analytical Results 

We now interpret the numerical results described above in 
terms of simple analytic formulae. Details may be found in 
the Appendix [X] Here we merely summarize the relevant 
results. 

In an axisymmetric force-free electromagnetic configu- 
ration, the drift Lorentz factor 7 can be described quite well 




Figure 6. Transverse variations of the field-line invariants, Q 
and |/|, for model A as a function of the magnetic field stream 
function P. In each panel four curves are shown: the star-disk 
surface (solid line), r = 10 2 (dotted), r = 10 4 (short-dashed), and 
r = 10 6 (long-dashed). Over a range of six orders of magnitude 
in distance, the values of Q and \I\ are conserved to better than 
12%. 



by the following analytic formula (see Appendix IA4. II) . 



-.9 -.2 "r 9 1 



V 71 
where 71 and 72 are given by 



7i 
72 



(SIR) 



12 



21 1/2 



fXR, 



1/2 



(13) 

(14) 
(15) 



where the last equality in (|14[) holds for Q.R 2> 1. Here, R — 
r sin 9 is the cylindrical radius, Q is the rotation frequency of 
the local field line, R c is the poloidal radius of curvature of 
the field line, and C is a numerical factor of order unity that 
depends on the field line rotation profile (see Appendix lA7|) . 



C: 



3 + 



d log Q. 
dlogR 



1/2 



(16) 



In the jet region we have Q = const and, therefore, we expect 
C « >/3 « 1.7. (17) 



As we will see in £13.31 in the simulations we find values of C 
slightly below this value because Q slightly decreases with 
increasing R toward the edge of the jet (due to numerical 
diffusion, see Fig. [6j) . 

Equation (|13p gives the drift speed of an infinitely 
magnetized magnetodynamic, or force-free, flow. One might 
worry that the velocity of a fluid carried along with such a 
flow will be very different. In Appendix|B]we show that any 
such fluid has only a slightly modified Lorentz factor relative 
to the drift speed, in the limit when the fluid is massless, 
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Figure 7. Similar to Fig. [3] but for model B, which corresponds to V = 1, Qo = 0.25, /3 = 1.5. Field lines in this model show faster 
collimation than in our fiducial model A. This is because the confining magnetic pressure of the disc wind falls off more slowly with 
distance - as r~ 2 here instead of as r~ 2,5 in model A. The difference in collimation is not due to hoop-stress or rotation. Note that the 
maximum Lorentz factor at any distance coincides with the jet-wind boundary, in contrast to what is seen in model A. 



i.e., a oo. Therefore, for all practical purposes we can 
assume that the fluid Lorentz factor is given by eq. (|13[1 . 

Since 7 2 is the harmonic sum of two terms, the value 
of 7 is determined by whichever of the two quantities, 71 
and 72, is smaller. Close to the central compact star, 71 
is smaller, and the first term in equation (|13p dominates. 
Thus, for a given rotation frequency Qo of the field lines in 
the jet (determined by the spin of the compact object), the 
Lorentz factor increases linearly with dis t ance from th e jet 
rotation axis dBlandford fc Znaiek|[l977l ; iBeskinl fl997l) . In 
this well-known regime, which we refer to as the first accel- 
eration regime, a faster compact object spin leads to faster 
acceleration along the jet. Also, for a given rotation, the out- 
ermost field lines in the jet, which emerge from the equator 



of the star, have the largest acceleration and largest 7 at 
any given z. 

The second term in equation (|13|l represents a slower 
regime of acceleration, which we refer to as the second ac- 
celeration regime. It is present only for certain field geome- 
tries and is generally realized only at large distances from 
th e compact obje c t. For the self-similar solutions described 
in lNaravan et all (|2007h . this acceleration regime is impor- 
tant only if the self-similar index v < 1. Since model A has 
v = 0.75, this term is important for our simulation. Note 
that models with v 1 are astrophysically the most inter- 
esting and relevant ( H5.ip . so it is important to understand 
the second acceleration regime. A feature of the second ac- 
celeration regime is that the Lorentz factor does not depend 
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explicitly on the field line rotation frequency, but is deter- 
mine d purely by the lo cal poloidal curvature of the field 
line (jBeskin et al] 12004 ) • 

Moreover, as we saw earlier, the 
poloidal structure of the field line is itself largely indepen- 
dent of rotation. 

Let us ignore the small distortion of the field line shape 
caused by rotation (Fig. [5} , and take the shape to be given 
by the initial nonrotating solution (see Appendix I A2[) : 



z oc R 



2/(2-") 



-1//2 



(18) 



The latter scaling, shown by the dashed line in Fig. [5] pro- 
vides a good description of the field line shape. Using this 
scaling, we can evaluate 71 and 72 in the jet using equa- 
tions HU) and (JT5j| (see Appendix lA4.2|) . 



71 » Q,o r sin6>, 

72 « k/0, 



(19) 
(20) 



where « = 1C I \J (2 — v)v does not have any explicit depen- 
dence on fi or position. This gives the following scaling of 
the Lorentz factor along field lines, 



71 oc r 

72 oc r 



l-v/2 

; 

v/2 



(21) 
(22) 



Close to the central star, 72 is always larger than 71, and 
thus the jet 7 is determined by 71. With increasing distance 
along a field line, 71 and 72 grow at different rates. If v > 1, 
72 rises more rapidly than 71 and the Lorentz factor of the 
jet is always determined only by 71 (e.g., model B below, 
which has v = 1). However, for v < 1 (model A and model 
C), 72 rises more slowly than 71 and takes over at a certain 
distance from the star. This corresponds to the second slower 
acceleration regime seen in Fig. [4] 

Figure [4] shows a comparison of the actual 7 measured 
in the model A simulation with the prediction from the an- 
alytic formula JT3]). We set C = V3 w 1.7 (see eq. ([TTJ) and 
the next subsection). We find that the analytic formula for 
7 agrees remarkably well with the numerical results. The 
formula gives the correct slopes and reproduces the distance 
at which the break between the two acceleration regimes 
occurs. 

As we see from Figure [4j the second regime of Lorentz 
factor growth is most prominent along field lines originat- 
ing closer to the equator of the compact object. This is the 
reason for the 'anomalous' development of a slower-moving 
sheath surrounding a faster-moving structured jet spine that 
we mentioned in £|3 . 1 1 See ij3.4l for more detail. 

The cause for the slight deviation of the magnetic pres- 
sure from the r -2 ' 5 power-law behaviour, as seen in Fig- 



ure[S] is discussed in Appendix lA6l We show that, for v < 1, 
the magnetic pressure shows a broken power-law behaviour 
along field lines, 



,2(v-a) 



r < 
r > 



(23) 



The break radius r tr is the same as the radius where the 
jet acceleration switches from the first regime (71) to the 
second (72). The power-law indices in (f2"3")l as well as the 
predicted break radius r tr ~ 7 x 10 3 are consistent with the 
results shown in Fig. [5] Note that the confining pressure of 
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Figure 8. Similar to Fig. [4] but for model B. The analytic fit 
for the Lorentz factor uses the same value of C = \/3 as before. 
Compared to model A, we see that model B lacks the second 
regime of acceleration. Correspondingly, there is no fast jet core 
present. 



the wind (along any field line originating in the disc suffi- 
ciently far from the compact object) follows a single power- 
law, r 2 '"" 2 ', at all distances from the compact object (see 

E3J. 



3.3 Dependence of the Results on Model Details 

In order to explore which features of the results described 
above are generic and which are particular to model A, we 
simulated a wide range of models with v varying from 0.5 to 
1.25. We find that model A is representative of most models 
with v < 1. In particular, all these models show the two 
regimes of Lorentz factor growth (|2ip and (|22[) . Similarly, 
we find that model B, which has v = 1 and is described 
below, is representative of models with v ^ 1. 

Model B has field lines with a parabolic shape, as we 
expect from equation (|18|l . Figures [7] [8] and [9] show results 
corresponding to this model. The jet acceleration is always 
in the first regime and the Lorentz factor of the jet is deter- 
mined only by 71 . Consequently, the maximum acceleration 
always occurs for the field line at the jet-wind boundary. 
This is obvious in Fig. [7] where we see that the maximum 
Lorentz factor coincides with the maximum in the enclosed 
current. Also, in Fig. [S] we see that the Lorentz factor of 
the field line with 9{ p = 20° is always smaller than that of 



the line with 8{ p 



64°. Our model B simulation is well- 



converged, with Q, and / preserved along the field lines to 
better than 15%, even though the simulation domain ex- 
tends over six orders of magnitude in distance. 

In model A, it was the presence of the second regime of 
Lorentz factor growth that was responsible for the develop- 
ment of a faster jet core. This regime is absent in model B 
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Figure 9. Similar to Fig. [5] but for model B. Note the excellent 
agreement between the numerical quantities and the analytical 
estimates. 
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Figure 10. The values of the factor C as determined from simula- 
tions of jets with different choices of the field geometry threading 
the star, the spin of the compact object, and the rotation pro- 
file of the disc. Crosses and diamonds correspond to simulations 
with Co = 0.25 and 0.1, respectively. The scatter is mostly due 
to variations of C from one field line to another (see text for 
more detail). For 0.5 u $J 0.8, the numerical results lie within 
the expected analytical range (125 b . between the solid and dotted 
horizontal lines. For v > 0.9 we are not able to reliably deter- 
mine the value of C from the simulation (an error bar is shown). 
However, for these values of u, the second regime of acceleration 
does not operate for any distance r < 10 6 , and the value of C is 
unimportant. Therefore the same value of C can be used here as 
well. 



(compare Fig. [S] and Fig. [SJ. Indeed, everything is much 
simpler in model B. For instance, Fig. [8] shows that the an- 
alytic formula for the Lorentz factor accurately reproduces 
the numerical profile, and Fig. [5] shows that both the field 
line shape and the comoving pressure accurately follow the 
predicted dependencies, 9 oc r~ 1//2 , p m ag oc r~ 2 . We obtain 
this kind of close agreement for all models with v ^ 1. 

We have investigated the sensitivity of the models to the 
rotation profile in the disc (the value of /3), the magnitude 



of the stellar spin (the value of firj), an d the geometry of 
the field threading the star. For v ranging from 0.5 to 1.25 
we tried different values of these parameters. In particular, 
we have done simulations with a uniform rotation velocity 
in the dis c, i.e., (3 = 1, which co rresponds to the self-similar 
model of iNaravan et all l|2007l ). and we have tried both a 
monopole field and a uniform vertical field threading the 
star. We find that these changes do not noticeably affect 
the jet; in particular, the field line shape changes negligibly. 
We have also investigated the effect of a slower stellar spin: 
flo = 0.1. We find that the field line shape stays very close 
to that of the nonrotating solution so long as v ^ 1, but 
changes logarithmically for v < 1, as in model A. 

We were particularly interested to see how well the gen- 
eral formula for the Lorentz factor (|13p performs for the 
range of models we considered. Since fi in the jet region is 
not perfectly constant due to inevitably present numerical 
diffusion (Fig. [S| , 



1<^<0, 



(24) 



d log R 

we expect a range of values for the factor C (eg. I16[) . 

V2 < C < VS. (25) 

The upper bound C = y/% (eq. \YJ\ is the analytical value 
for the case fl = const. Figure [101 shows that for 0.5 ^ v ^ 
0.8 the best-fit values of the factor C for various field lines 
threading the star are within the expected range (125[) . for 
all models we considered. For v > 0.9 the second regime of 
the Lorentz factor (| 15p is unimportant (it is realized if at all 
only at greater distances than are of interest to us), and so 
the value of C is irrelevant. Thus, we can use the analytical 
value of C (|17p with equation (I13|) for all values of v in the 
range 0.5 ^ v ^ 1.25, for all physically relevant values of 
flo, and for any value of /3 between 1 and 1.5 (we have not 
explored other values). In all cases, for most of the jet (for 
field lines with ^ 9{ p ^ 80°), fl is constant to within a few 
percent and the Lorentz factor predicted by equation (|13|) 
agrees with the numerical results to better than 10% (see, 
e.g., Figs. |4]and [8j). 

3.4 Collimation and Transverse Structure of Jets 

Figures 3] EI [HI and OH show the behaviour of various quan- 
tities along field lines in models A and B. We now consider 
how these quantities vary across the jet at a given distance 
from the compact object. The results are shown in Fig. 1111 
The upper panel of Fig. [11] shows the angular profile 
of the Lorentz factor at various distances from the central 
object for each of the three models, A, B, C. Consider the 
simplest of the three models, model B, which has v = 1. As 
Fig. [8] shows, in this case the Lorentz factor is simply equal 
to Q.R at large distances. Since the field lines in the 'jet' 
region of the outflow are all connected to the rigidly rotating 
compact object at the center (see Fig. [2]), all of them have 
the same Q = flo- Therefore, we expect 7 = QoR ~ Q.ar6. 
This linear increase of 7 with 9 at a fixed r terminates at 
the edge of the jet, 6 = 9j. At the jet boundary we have 

9 = 9^Jl, 7 « 2^ = M « 0.35V?. (26) 
V r 9 j 9 j 

Beyond the edge of the jet, the field lines are attached to 
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Figure 11. Upper Panel: Angular dependence of the Lorentz fac- 
tor 7 for the fiducial model A (y = 0.75, solid lines), model B (y = 
1, dashed lines), and model C {y = 0.6, dot-dashed lines). For each 
model, three distances from the compact object are shown. From 
bottom to top, they are r = 1.25 X 10 5 , 5 X 10 5 , 2 X 10 6 . For a 
maximally spinning BH of mass M = 3Mq, these distances corre- 
spond to r = 6 X 10 10 , 2 X 10 11 and 9 X 10 11 cm; for a neutron star 
of radius 10 km, the distances in centimeters are approximately 
half these values. Asterisks indicate the boundary between the 
jet and the wind. Note that for all models and distances shown, 
7 is ~ 100 or larger. Note also that, for models A and C, the 
maximum value of 7 occurs inside the jet, whereas for model B, 
it occurs at the boundary between the jet and the wind. Lower 
Panel: Angular dependence of the jet energy flux per unit solid 
angle, dP/dui, for the same models and distances. The energy flux 
is maximum at the jet-wind boundary for all three models. 



the disc, where £1 falls rapidly with increasing radius. There- 
fore, the Lorentz factor decreases quickly. Appendix IA4.2I 
shows that the expected dependence is 7 oc 9~ 2 . The dashed 
lines in Fig. [11] confirm the scalings of 7 with 9 both inside 
and outside the jet. They also show how 9j and 7 vary with 
distance from the compact object. 

Consider next model A, with v — 0.75. Now we expect 
the jet angle to scale as 



V2_ 

W 2 



r 3/8 ' 



(27) 



This is approximately verified in Figs. [S] and 1111 However, 
the agreement is not perfect because field lines open out 
slightly at large distances relative to the analytical approx- 
imation of the field line shape. 

The Lorentz factor profile has a more complicated be- 
haviour in this model because of the presence of two different 
regimes of acceleration (Fig. Q. Close to the axis, the field 
lines are in the first acceleration regime, where 7 ~ £loR and 
the behaviour is the same as in model B, i.e. 7 oc 9. How- 
ever, at an angle 9 = 9 m ('m' for maximum), we begin to 
see field lines that have switched to the second acceleration 
regime, and beyond this point the Lorentz factor decreases 
with increasing 9. Thus we have 



7: 



n r<9, 
3.8/0, 



< e < e m , 

dm < 9 < 0j- 



(28) 



The coefficient 3.8 is obtained from the simulations, but it is 
close to the analytical value of k = 3.6. (The small difference 
is because of the slight opening up of the field lines in the 
second acceleration regime). The angle corresponding to the 
maximum Lorentz factor is 



3.8 
tt r 



3£ 



(29) 



and the value of the maximum Lorentz factor is (eq. 1 13 p 



7" 



—=Q, a r9 rl 
V2 



%/TMhr « 0.7Vr. 



(30) 



Interestingly, over the entire range of angles from 9 m to 
9j, we have the simple scaling 

3.8 



7 : 



< 9 < . 



(31) 



Note also that the coefficient in this relation is larger than 
the one corresponding to model B (eq. I26[l . Thus, for the 
same Lorentz factor, the jet in model A has at least six 
times larger opening angle than the jet in model B. Beyond 
the edge of the jet, in the wind region, the Lorentz factor 
drops rapidly as 9~ 2 5 . 

We can now make a general prediction for how the peak 
Lorentz factor j m scales with distance from the compact 
object. For a maximally spinning BH, 7 m attained at a dis- 
tance r is likely to lie in the range bounded by models A 
and B feqs. [Poland 1317)) : 



7m (r) = (0.3-0.7)^- 



(32) 



Note that this formula gives the maximum value of the 
Lorentz factor over all field lines, as opposed to a value of 
the Lorentz factor along a single field line. 

As we see from Fig. 1111 model C {v = 0.6) is a more 
extreme version of model A. The jet is significantly wider 



© 2008 RAS, MNRAS OOO.ITH251 



Simulations of Ultrarelativistic Magnetodynamic GRB Jets 13 



(8j ~ 0.2 radians), the maximum Lorentz factor is much 
larger (~ 5000 at r = 2 x 10 6 ), and the maximum of 7 oc- 
curs at 6 <C 6j . Unfortunately, our numerical results for this 
model show large deviations from the analytical model de- 
scribed in the Appendix. In particular, the poloidal field line 
shape is nearly monopolar at large distance, i.e., 6 is nearly 
independent of r along each field line, instead of behaving as 
~ r -0 ' 3 (see eq.QHUlE] 

Qualitatively, however, model C is 
similar to model A. Model C is well-converged with £1 and / 
preserved along the field lines to better than 15% through- 
out the simulation domain. 

We consider next the power output of the jet. We define 
the angular density of electromagnetic energy flux as 



(33) 



where S is the Poynting vector and dui = dtp d8 sin 6 is the 
solid angle. As we show in Appendix lA8l and verify in Fig. 1111 
the power output grows quadratically with distance from the 
jet axis for all models, 



dP 

dui 



1 2 2v n 2 



(34) 



reaches its maximum at the jet boundary, indicated by as- 
terisks in Fig. 1111 



dP 

did 



0.02 



2-4/3/iz 



(35) 



Equa- 



and falls off rapidly in the wind region as 6 
tion p5p illustrates that jets with a smaller opening angle 
have a larger power output per unit solid angle. The steep- 
est decline of angular power in the wind region is 8~ 10 in 
model C, followed by less steep declines of 6~ 6 in model A 
and 8~ 4 in model B. This behaviour can be seen in the 
lower panel of Fig. 1111 The angular power output profile in 
model C does not evolve with distance since the opening 
angle of the jet is nearly independent of distance. 

The total power output in the jet and the wind is (see 
Appendix [ 



2 



32 



i3z 



208/1/ - 1) 32(/3/^-l) : 



(36) 
(37) 



where B r is the radial magnetic field strength near the BH. 
The total power is independent of the distance at which it is 
evaluated, which is a manifestation of the fact that energy 
flows along poloidal field lines. The total power output in 
the jet is the same for all models, and the total power in the 
wind varies from model to model. The most energetic wind 
is in model B and carries twice as much power as the jet. 
In model A, the wind and the jet have equal power outputs, 
and in model C, the power output in the wind is two thirds 
of the power in the jet. 

To obtain the total power of the jet in physical (cgs) 
units, we specialize to an astrophysical system with a BH of 



1 We note that the monopole has a low acceleration effi- 
ciency for a finite value of m agnetization a llBeskin et al,|[T998l ; 
iBogovalov fc Tsineano"slll999l) . Therefore, for a finite magnetiza- 
tion, the acceleration efficiency might be also low for model C. 



mass M, radial field strength near the BH B r , and angular 
rotation frequency of field lines ilo- We obtain 



P 1CZ = ^«oS>oC 



1.8 x 10 D 



go 



M 



\ 10 15 G J V 3M P 



(38) 



where ro w r g = GM/c? for a rapidly spinning BH. For 
~ fimax = 0.25c/r 9 , B r ~ 1Q 15 G, 



characteristic values, fi 
M ~ SMq, and taking the typical duration of a long GRB 
~ 10 — 100 s, the model predicts a total energy output of 
10 51 — 10 52 erg, whic h is compara ble to the energy output 
inferred for GRB jets (Piran 2005). We note that the actual 
value of the magnetic field B r might be higher since the ob- 
servations only account for a fraction of the electromagnetic 
energy flux. Recent GRMHD simulations suggest a value of 
10 le G (|McKinnevll2005bD . 



4 COMPARISON TO OTHER WORK 

Since the observed ener gy output of long GRBs is 10 51 — 
10 52 erg (|Piranl l2005h . any model of long GRBs re- 
quires a central engine capable of supplying this copious 
amount of energy. In the absence of magnetic fields, a 
possible energy so urce is annihilation of neutrinos from 
the accretion disc (iKohri et a.1.1 120051 : IChen fc Beloborodovl 
120071 ; iKawanaka fc Mineshigej I2007T ). Attempts to include 
the neutrino physics self-consistently ha ve so far not 
succe e ded in produc i ng re lativistic jets (Nagatak i et al.l 
120071 ; iTakiwaki et al.l I2007T I. However, an ad-hoc quasi- 
isotropic energy deposition of 10 50 — 10 52 ergs into the 
polar regions seems to lead t o jets with Loren t z fac- 
tors ~ 100 jlAlov et all l200fj|; IZhang et all l2003t |2004| ; 



iMorsonv et alJl2007l : IWang et al.ll2007h . The jets so formed 
are found to be stable in two and three dimensions, to accel- 
erate via the conversion of internal energy into bulk kinetic 
energy, and to collimate to angles < 5° as a result of inter- 
acting with the dense stellar envelope. We note that similar 
hydrodynamic simulations of short GRB jets have attained 
much higher Lorentz factors ~ 700 due to the absence of 
a ste llar envelope that the jet has to penetrate l|Alov et al.l 
l2005h . 

The inclusion of magnetic fields enables the 
extraction of energy from spinning BHs via the 
Blandford-Znajek mechanism (|Blandford fc Znaiekl Il977l ; 
iKomissarov fc McKinnevI l2007h and f rom accretion discs 
via t he action of magnetic torques (Blandfo rd fc Pavnel 
1982). The paradigm of electromagnetically powered jets 
is particularly appealing since it is able to self-consistently 
repr oduce the observe d jet energetics of long GRBs (see i]3.4l 
and lMcKinnevl[2005bl ). without any need for ad-hoc energy 
deposition. Numerical simulations, within the framework of 
general relativistic magnetohydrodynamics (GRMHD), of 
magnetized accretion discs aro und spinning BHs naturally 
produce mildly re l ativistic jets dMcKinnev fc Gammieil2004l; 
McKinnevI l2005bl; lDe"Vilriers et al.ll2005l; lHawlev fc Krolikl 



20061 ; iBeckwith " et al.ll2007l ; iBarkov fc Komissarovl 120071) as 
well as relativistic jets with Lorentz factors ~ 10 ( McKinnevI 
l2006bh . 

Our results are consistent with those of iMcKinnevI 
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(2006b) of a highly magnetized jet that is supported by 
the corona and disc wind up to a radius r ~ few x 10 2 
by which the Lorentz factor is ~ 10 similar to expected by 
our v — 3/4 model. They did not include a dense stellar 
envelope. Beyond the distance r ~ few x 10 2 the disk wind 
no longer confines the jet, which proceeds to open up and 
become conical as it passes the fast magnetosonic surface. 
Beyond the fast magnetosonic surface they find the jet is 
no longer efficiently accelerated, as consistent with expec- 
tations of unconfined, conical MHD solutions dBeskin et al.l 
19981; iBogoyalov fc Tsinganosl Il999h - The jet simulated by 
McKinnevi (12006b) shows a mild hollow core in 7 but shows 



no hollow core in energy flux, although this may be a result 
of numerical diffusion or significant time-dependence within 
the jet. 

iBucciantini et all j|2006l . I2007I . l200St ); 

iKomissarov fc Barkovl (|2007 ) studied the formation 
and propagation of relativistic jets from neutron stars as a 
model for core-collapse-driven GRB jets. They were unable 
to study the formation and propagation of ultrarelativistic 
jets likely due to computational constraints and their model 
setup. While our model is somewhat idealized compared 
to those studies, we demonstrate the generic process of 
magnetic-driven acceleration to ultrarelativistic speeds and 
thu s extend their simulatio n s of ma gnetar-driven GRB jets. 

iMcKinnev fc Naravanl (|2007al lbh found that the jets 
in their simulations are collimated by the pressure sup- 
port from a surrounding ambient medium, such as a disc 
corona/wind or stellar envelope, rather than being self- 
collimated. Thus the treatment of the medium that con - 
fines the jet appears to be cr ucial. IKomissarov et (|2007h : 
iBarkov fc Komissarovl l|2008f l modelled the action of such an 
ambient medium by introducing a rigid wall at the outer 
boundary of the jet and studied the magnetic accelera- 
tion. They obtained solutions with Lorentz factors of ~ 300 
with an efficient conversion of magnetic to kinetic energy. 
In the current work, instead of keeping the shape of the 
jet boundary fixed via a rigid wall, we prescribe the ambi- 
ent pressure profile, so that the shape of the jet boundary 
can self-adjust in response to pr essure changes inside the 
jet (|McKinnev fc Naravanl 1200781 ). We believe that this is 
a more natural boundary condition for modeling GRB jets 
confined by the pressure of the stellar envelope. 

The magnetodynamical, or force-free, approximation 
may provide a good approximation to the field geometry 
even in the mass-loaded MHD regime as long as the flow is 
far from the singular monopole case of an unconfined flow. 
iFendt fc Quvedl ( 20041 ) used this fact to study ultrarelativis- 
tic jets in the GRB context. As in our work, they determine 
the field line geometry in the magnetodynamical regime, but 
they use energy conservation for the particles to determine 
an approximate MHD solution and approximate efficiency 
of conversion from magnetic to kinetic energy. 

The conversion efficiency for fully self-consistent highly 
relativistic MHD solutions has been studied only for a lim- 
ited number of field geometries. High conversion efficiency 
was found for a parabolic v = 1 solution and an intermedi- 
ate v — 2/3 solut ion, but not for th e singular monopole 

v = solution dBeskin et aL 19981 ; iBeskin fc Nokhrinal 

120061 ; IBarkov fc Komissarovl 120081 ). Following this work 
we plan to include particle rest-mass and systematically 
study the efficiency of particle acceleration in highly rel- 



ativistic self-consistent MHD solutions. For this we will 
use the same numerical scheme as in the present nu- 
merical work but optimized f or the ultrarelativistic MHD 
regime llGammie et ail 120031 ; iMignone fc McKinnevi 120071 ; 
iTchekhovskov et alj|2007l ). 

Since we assume axial symmetry in this study, our simu- 
lations cannot address the question of jet stability to the 3D 
kink (screw) mode. Unlike the hydrodynamic case, we are 
not aware of any studies of relativistically magnetized jets 
in 3D that address this issue. However, all jets in this paper 
have B v « —QRB r (see eq, IA42ll . Thus, they marginal ly sat- 
isfy the stability criterion of iTomimatsu et alj 0001), sug- 
gesting that our jets are marginally stable to kink instabil- 
ity. Also, the spontaneous development of the spine-sheath 
structure in models with v < 1 may be naturally stab ilizing 
|Mizuno et al.ll2007l ; lHardedl2007l ; Hardee et al1l2007l ). 



5 ASTROPHYSICAL APPLICATION 

Since we consider magnetodynamic, or force-free, jets in 
this paper, we are not able to study the effects of mass- 
loading of the jet. However, we expect the main properties 
of our infinitely magnetized jets to carry over to mass-loaded 
jets, provided the latter are electromagnetically dominated. 
Mass-loaded jets stay electromagnetically dominated as long 
as the Lor entz factor is well below th e initial magnetization a 
of the jet l|Beskin fc Nokhrinall2006l ) , where a is the ratio of 
electromagnetic energy flux to mass energy flux at the base 
of the jet. We assume that this condition is satisfied and pro- 
ceed to apply our results to GRBs and other astrophysical 
systems with relativistic jets. 



5.1 Application to Long GRBs 

The first question of interest is what sets the terminal 
Lorentz factor of a relativistic jet. We have shown in this 
paper that 7 increases with distance from the compact ob- 
ject as 7 ~ l/0j ~ r" /2 . The value of v is determined by 
the radial dependence of the confining pressure: if pressure 
varies as r~ a ~ r~ 2 ' 5 , then v = 2 - (a/2) ~ 0.75. 

In the context of the collapsar or magnetar model 
of GRBs, the confining pressure is primarily due to the 
stellar envelope, and hence acceleration is expected to 
continue only until the jet leaves the star. Once out- 
side, the jet loses pressure support and the magnetic field 
configuration will probably become conical (monopolar). 
This geometry is inefficient for accele rating particles to 
Lorentz factors larger than 7m ax ~ a ' l|Beskin et al. I ll998l ; 
iBogovalov fc Tsinganosl Il999h . We note that mildly rela- 
tivistic edges of the jet will expand quasi-spherically after 
the jet loses pressure support, while the ultrarelativistic jet 
core will open up at most conically l|Lvutikov fc Blandfordl 
2003). In fact, if the jet loses pressure support at a distance 
r ~ r a (which corresponds to the radius of the progenitor 
star), the opening angle of the jet core will stay approxi- 
mately constant until the jet reaches a much larger distance 



> 10 r s (for our models A and C). We note 



also that current-driven instabilities may set in when the 
jet loses pressure support, and much of the elect romagnetic 
energy may be conve rted into thermal energy (|McKinnevl 
2006b; L vutikovl 120061 ). These additional topics are beyond 
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the scope of the present paper and require simulations that 
include a loss of pressure support at large radius to model 
the stellar surface. 

From the above discussion, we expect that the termi- 
nal Lorentz factor of the jet in a long GRB is determined 
primarily by the size of the progenitor star. Figure HTl shows 
that, for a BH of a few solar masses and stellar radii in 
the range few x 10 10 cm to 10 12 cm, the Lorentz factor of 
the jet is expected to be between about 100 and 5000. Such 
masses and stellar radii are just what we expect for GRB 
progenitors, and the calculated Lorentz factors are perfectly 
consistent with the values of 7 inferred from GRB observa- 
tions. 

For our fiducial model, the size of the star determines 
the terminal Lorentz factor as long as the initial magneti- 
zation of the jet is sufficiently high, a > 10 3 , but not infi- 
nite, a < 10 9 . The first condition ensures that the jet inside 
the star is well-described by the force-free approximation. 
However, once the jet breaks out of the star and its geom- 
etry becomes monopolar, the effects of finite magnetization 
will set in, provided the second condition is satisfied. Ac- 
celeration in such monopolar field is ineffective, therefore 
the growth of the Lorentz factor stops. Estimates of baryon 
loading of magnetized GRB jets from black hole accretion 
systems are somewhat uncertain. We can estimate the initial 
magnetization of a GRB jet as the ratio of electromagnetic 
jet powe r output ([38]) to the rate o f baryon mass- load ing of 
the jet ()Levinson fc Eichler|[2003al : iMcKinnevI l2005al . their 
cq. A10), which gives a value of initial jet magnetization 
a ~ 10 3 for characteristic values of the accretion system 
parameters. However, given the uncertainty in these param- 
eters, the actual value of magnetization could be an order 
of magnitude lower or higher. Note that magnetically accel- 
erated jets behave very differently from relativistic fireballs. 
Beyond the star fireballs adiabatically expand leading to 
7 oc r until nearly all thermal energy is exhausted, and so 
there is no spatial scale that sets a terminal Lorentz factor. 
Fireball models require the baryon-loading of the jet be fine- 
tuned to obtain the observed Lorentz factor and opening 
angle. However, in the MHD case, acceleration essentially 
ceases at the stellar surface for a wide range of values of a. 

We discuss next the degree of collimation of the jet. 
The three models we have described in this paper form a 
one-parameter sequence. Model B has the most highly col- 
limated jet and model C has the least collimated jet, while 
our fiducial model A lies in-between. From Figure [11] we 
see that the models produce jet opening angles in the range 
10 -3 — 0.2 radians, which is consistent with the ty pical open- 
ing angles observed in long GRBs ~ 0.03—1 rad (Frai l et al.l 
l200ll ; IZeh et~ai1l2006l ). In making this comparison, we are as- 
suming that the opening angle of the jet is set by its value 
when the jet emerges from the star. 

Observations of GRB afterglows often indicate an 
achromatic 'jet break' in the light curve a day or two after 
the initial burst. These breaks suggest that the opening an- 
gle 8j of the jet is substantially larger than the initial beam- 
ing angle I/7 of the GRB. Roughly, it seems that GRBs have 
y0j > 10. In this regard, the three models described in this 
paper have very different properties. Model B (y — 1) has 
"f6j < 1. This model would predict a jet break immediately 
at the end of the prompt GRB, which is very different from 
what is seen. Model A (u — 0.75) has jOj ~ 4 for the plasma 



at the boundary of the jet. This model has larger values of 
7 in the interior, so it predicts a range of values of jOj from 
4 up to about 15. Model C (y = 0.6) is even more extreme, 
giving values of j8j ~ 25 — 950. These models span a range 
that is more than wide enough to match observations. Our 
fiducial model with v — 0.75 appears to be most consistent 
with the data. 

Note that the different models correspond to different 
radial profiles of pressure in the confining medium: model B 
has pressure varying as r~ 2 , model A goes as r~ 2 ' 5 , and 
model C as r~ 2 ' 8 . We selected model A as our fiducial model 
because a superficial study of GRB progenitor star models 
suggested that the pressure probably varies as r~ 2 ' 5 . How- 
ever, in view of the fact that jet properties depend strongly 
on the value of the index, a more detailed study of this point 
is warranted. 

Finally, we consider the jet power. The total energy 
output in a GRB jet is given by equation (|38p . This as- 
sumes that the power delivered by the disc wind is de- 
posited entirely into the stellar envelope and so does not 
contribute to the jet, although the mass-loading of field 
lines near the jet-wind boundary proba bly depends upon 
the amount of turbulence in the interface l|McKinnevll2006bl ; 
iMorsonv et al.ll2007l : IWan g et ai1l2007f) and non-ideal diffu- 
sion physics (|Levinson fc Eich lcr 2003bl; IMcKinnevI [2005a; 
iRossi et al.ll2006h . For characteristic values of the magnetic 
field strength B r near the compact object and the angu- 
lar frequency Oo of the compact object, the energy released 
over the typical duration of a long GRB is about 10 5 



10 J 



erg, which is consistent with the observed power (|Piranl 
2005). Other things being equal, the e nergy scales as Qq. 
This scaling (jBlandford fc Znaiek|[l977h would certainly ap- 
ply to a magnetar. In the case of an accreting black hole, 
the magnetic field strength near the black hole may itself 
increase rapidly with Qo (for a fixed mass accretion rate), 
and this may give a steeper dependence of jet power on fi 
|McKinnevll2005bh . 

Figure [TT] indicates that the electromagnetic energy flux 
in the jet has a substantial variation across the jet. The 
power is very low near the jet axis, and most of the energy 
comes out near the jet boundary. Note in particular that, 
for the observationally relevant model A, the maximum jet 
power does not coincide with the maximum Lorentz fac- 
tor (see Fig. [12}. This unusual behaviour is the result of 
the second acceleration regime which we discussed in H3.4I 
(see Figs. [3] [4] 11 1 pi . Such a ring-like shaped distribution of 
Lorentz factor in a fireball naturally leads to the Amati re- 
lation of an observed correlation between de-redshifted peak 
frequency in the GR B spectrum and the isotrop ic equivalent 
luminosity of GRBs (|Eichler fc Levinsonll2004r ). Regardless, 
the fact that the jet power comes out along a hollow cone, 
and not as a uniform ly fil led cone as assumed for example by 
iRhoadsl (|l997l , ll999l ) and lModerski et al.l < |2000h . is likely to 
have observational consequences for both the prom pt emis- 
sion from a GRB and its afterglow (|Granoj|2005l ). We are 
assuming, of course, that the electromagnetic power which 
we calculate from our model is directly proportional to the 
obs erved radiative pow er (both prompt and afterglow). 

IRossi et all {2002) suggested the interesting possibility 
that GRB jets may be 'structured,' with the energy flux per 
unit solid angle, dP/dui, having a flat core and the power 
falling off as ~ 8~ 2 outside the core. We have already seen 
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Figure 12. Face-on view of the lateral structure of the jet in model A at a distance of r = 5 X 10 5 (2 X 10 11 cm for a maximally spinning 
BH of mass M = 3Mq). The left panel shows the energy flux of the jet per unit solid angle, dP/dw, and the right panel shows the 
Lorentz factor, 7. Red colour indicates large and blue colour small values of the respective quantities. The dashed line shows the jet- wind 
boundary. Note that the maximum value of the energy flux, dP/dw, occurs at the jet- wind boundary, while 7 is maximal inside the jet. 



that our electromagnetic jets do not have a fiat core. In ad- 
dition, we find that the power outside the jet, in the wind re- 
gion, falls off very steeply. The variation is 6T 4 in the mildest 
case (model B) and is as steep as 8~ 10 in the most extreme 
case (model C). In addition, the disc wind in our idealized 
model is merely a proxy for a gaseous confining medium, 
which means that the electromagnetic power in this region 
may be even less than we estimate. Entrainment and insta- 
bilities within the jet may lead to a less sharp distribution 
at large angles. 



5.2 Application to AGN, X-ray Binaries, and 
Short GRBs 

In the case of accreting black holes in active galactic nu- 
clei (AGN), accreting neutron stars and black holes in X- 
ray binaries, or accreting black holes for short GRBs (which 
presumably form as a result of a coalescence of c ompact 
object binary systems, IPiranl 120051 ; iMeszarosI l2006h . there 
is no stellar envelope to confine the jet. Therefore, the 
only confining medium available is the wind from the in- 
ner regions of the accretion disc. The strength of these 
winds is known to be a strong function of the radia- 
tive efficiency of the disc. A radiatively inefficient disc 
(i.e., an advection-dominat ed accretion flow, AD AF), will 
have a strong disc wind (|Naravan fc Yil [l994), whereas 

a radiatively efficie nt disc (i.e., a standard thin disc, 

IShakura fc Sunvaevl 1 19731 ). will normally have a much 
weaker wind. Thus, jet acceleration should be more effec- 
tive in systems with ADAFs, viz., low- luminosity AGN and 
black hole binaries (BHBs) in the 'hard state' and the 'quies- 



cent state' |Naravan fe McClintockll2008l ). Indeed, observa- 
tions indicate that these systems are invariably radio-loud, 
whereas higher luminosity AGN and BHBs in the 'thermal 
state', which are powered by thin accretion discs, are of- 
ten radio quiet. On the other hand, some of the most ener- 
getic radio quasars are associated with high luminosity AGN 
which presumably have thin discs. It is unknown how such 
jets are confined, but large-scale force-free magnetic fields 
could replace the material support of the disk or wind. 

The terminal Lorentz factor of the jet in an ADAF or a 
short GRB system will depend on the distance out to which 
the disc wind is able to provide significant pressure support. 
Numerical GRMHD simulations of AD AF-like tori of ex - 
tent r ~ 40 around rapidly spinning BHs (McKinney 2006b) 
show that the disc wind is effective to a radius r ~ few x 10 2 , 
giving a terminal Lorentz factor 7 ~ few — 10 as long as 
a < 10 s . Such L orentz factors are c onsisten t with values 
inferr ed for AGN (Ijorstad et alj|2005h. B HBs (jFender et all 
l2004h . and short GRB jets (|Nakarl l2007l V Under some cir- 
cumstances neutr on star X-ray bina ries may produce jets 
similar to BHBs l|Fender et all l2004h . Short GRBs can be 
mass-loaded by neutron deposition and are limited to a < 
10 3 just as required. 

Once the confining effect of the disc wind ceases at 
r ~ 100, further acceleration will require some other confin- 
ing medium, e.g., an external interstellar medium, to take 
over. This would appear to be unlikely for a BHB or a short 
GRB, but it may work for some AGN. Let us assume that the 
jet stops collimating (becomes monopolar) beyond r ~ 100. 
The opening angles that we expect then lie in the range 
bounded by the most collimated model B and the least col- 
limated model C, 9j ~ 0.14 — 0.5 (eg. I27[l. i.e. substantially 
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larger than for long GRBs. These values are perfectly con- 
sistent with the jet opening angles inferred for the systems 
considered in this s ubsection (see the above references and 
IWatson et afll2006h . 



6 CONCLUSIONS 

By performing axisymmetric time-dependent numerical sim- 
ulations, we have studied highly magnetized, ultrarelativis- 
tic, magnetically-driven jets. In the context of the collapsar 
scenario of long GRBs, we obtain global stationary solu- 
tions of magnetodynamic, or force-free, jets confined by an 
external medium with a radial pressure profile motivated by 
models of GRB progenitor stars (see J2|. 

We find that both the size of the progenitor star and 
the radial variation of pressure in the envelope of the star 
determine the terminal Lorentz factor and the opening angle 
of the jet for a wide range of initial magnetization of the jet. 
At the radius where the jet breaks out of the star, i.e. at a 
distance r ~ 10 10 - 10 12 cm, the jets in our models attain 
bulk Lorentz factors 7 ~ 100 — 5000. These are currently the 
largest 7 attained in a numerical MHD jet simulation of long 
duration GRBs. The Lorentz factors we obtain are perfectly 
consistent with the values of 7 inferred from observations. 
The simulated jets have opening angles in the range 9j ~ 
10" 3 - 0.2 radians, in agreement with the typical opening 
angles observed in long duratio n GRBs (~ 0.03 — 1 rad, 
iFrail et al.ll200ll ; IZeh et al.ll2006h . 

For a maximally rotating black hole ora~ 1-ms magne- 
tar, a characteristic magnetic field near the compact object 
of 10 G, and a burst duration of 10 — 100 seconds, our 
simulated jets provide an energy output of 10 51 — 10 52 erg, 
comp arable to the power output inferred for GRB jets (|Piranl 
120051 ). 

The angular structure of the simulated jets is not uni- 
form. The jet power comes out in a hollow cone, peaking 
at the jet boundary (Fig. I12p . However, the Lorentz fac- 
tor peaks neither at the jet axis nor at the jet boundary 
but rather in between (Fig. I12|l . This nonuniform lateral jet 
structure may have observational consequences for both the 
prompt emission from GRBs and afterglows. 

To fully interpret simulation results, we derive a simple 
approximate analytical model (Appendix[S} which gives the 
scaling of 7 and 0j as a function of distance away from the 
compact object, as well as the variation of 7 and energy flux 
across the face of the jet. With these scalings we are able 
to understand all our simulation results, and we can predict 
the properties of highly magnetized jets in more general sit- 
uations. In particular, in most of our models we find that 
the maximum Lorentz factor at a distance z away from the 
compact object is j m ~ (z/ro) 1//2 , where ro is the central 
black hole/magnetar radius feq. I32|l . 

While the magnetodynamical regime that we have stud- 
ied allows us to establish the ability of highly magnetized 
flows to accelerate and collimate into ultrarelativistic jets, 
this approximation cannot be used to establish the efficiency 
of conversion from magnetic to kinetic energy. Following 
this work we plan to include particle rest-mass to study 
the properties of jets in the MHD regime and to determine 
the efficiency with which plasma is accelerated in the ul- 
trarelativistic MHD regime. Future simulations will focus 



on the time-dependent formation of ultrarelativistic magne- 
tized jets from the central engine and propa gation of the jet 
through a realistic stellar env elope (see, e.g.. lTakiwaki et al.l 
120071 ; iBucciantini et al]|2008t ) . Our present magnetodynam- 
ical and future cold MHD results should be a useful theoret- 
ical guide for understanding these more realistic and com- 
plicated simulations. 
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APPENDIX A: APPROXIMATE ANALYTICAL 
DESCRIPTION OF MAGNETODYNAMIC JETS 

Al Approximate Solution for the Flux Function 

An axisymmetric steady non-rotating force-free configura- 
tion is described in spherical co ordinates by a magne tic flux 
function P{r, 8) which satisfies (|Naravan et al.| [2007) 



1 8P 



2 d'P . 8_ 
T dr 2 m 88 \ sin6» 88 



= 0. 



(Al) 



Specializing to a self-similar configuration, we look for a so- 
lution of the form 

P(r,9)=r»p(6), (A2) 

where we limit ourselves to ^ v ^ 2 (|Naravan et al.ll2007t ) . 

Substituting (|A2I) into (|A1[I . we obtain the following 
differential equation for p(8), 

p"{8) - cot Op {8) + v(y - l)p(8) = 0, (A3) 

with boundary conditions, p(0) = and p(ji/2) = 1. The 
general solution is 

p{8) = ci/i(0) - c 2 cos 8 f 2 {8), (A4) 

where /i and /2 are hypergeometric functions, 



M") = 2l\[ 



(A5) 
(A6) 



and the constants ci and C2 are determined from the bound- 
ary conditions, 



Ci = 1, c 2 = 



vV (3/2 - v/2)T(v/2) 
r(l-i//2)r(!//2 + l/2)' 



(A7) 



For the particular cases of v — 1 and v = 0, we find 
that fi(8) = fzip) = 1, and the solution is very simple: 
p(8) = 1 — cos 8. Using this result as a guide, we write down 
the following approximate solution for the general case: 



P(r,8) w r"(l - cos( 



(A8) 



This turns out to be a good approximation to the exact 
solution for ^ v ^ 1.25. The relative error in P is less 
than 10% over this entire range of v. 

The most attractive feature of (|A8|) is its simplicity. Fur- 
thermore, although we began by considering a non-rotating 



configuration, the same solution turns out to be a reason- 
able approximation even for the rotating case, at least for 
rotations up to the maximum we consider, f2 ma x = 0.25. 
Thus, using (|A8|I , we can obtain a variety of useful, though 
approximate, results to describe force-free jets. 



A2 Field Line Shape 

By definition, the flux function P(r, 8) is constant along a 
field line. Therefore, it is straightforward to obtain from 
equation (|A8|) the shape of a field line as a function of the 
position of its footpoint, r = rf p , 8 = 8{ p : 



I. — cos 8 — ( — J (1 — cos 0f p ), 



which can be rewritten as 



sin - — 

2 



rfp 



-v/2 



(A9) 



(A10) 



We are primarily interested in the properties of the jet at 
large distances from the compact object, where field lines 
collimate and <C 1. In this limit, the field line shape is 
given by 

— 2sin^, (All) 

r fp ) 2 



or, inverting the formula, 



r 

r (p 



-2/v 



2 sin — - 



2/v 



In cylindrical coordinates, the field line shape is 



_R 

rip 

z 



l-v/2 



2 sin ■ 



2 ' 



2/(2-!-) 



2sin- 



-2/(2-") 



(A12) 

(A13) 
(A14) 



In the models we analyzed in the main text of the paper, 
we defined the 'jet' to include all field lines whose footpoints 
are on the compact central object, rf p = 1, ^ 6>f p < n/2, 
and the 'wind' to consist of the remaining field lines with 
footpoints on the disc, rf p > l,6>f p = n/2. Thus, the field 
line that emerges at rt p = 1, 8{ p — n/2 defines the boundary 
between the jet and the wind. The angle 9j of this field line 
at a distance r is 



0j w V2r 



-v/2 



(A15) 



This relation shows how collimation proceeds as a function 
of distance for a given value of v. Since the models of interest 
to us have 0.5 < v < 1, the collimation is gradual, going 
more slowly than 1/ 1 \fr. 

Field lines with 8 < 8j connect back to the central ob- 
ject (rfp = 1), and their foot-points are located at 

1 



r fp 



7fp 



V2> 



«S 8 < 



(A16) 



Field lines at larger angles connect back to the disc (#f p 
7r/2), and their foot-point radii are given by 



rfp 



2/v 



7fp 



tt/2, 



(A17) 
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A3 Poloidal Magnetic Field 



B r = r , 

B e = -ur v ' 2 tan{9/2). 



Given the magnetic flux function (1A8[) . we can determine the 
poloidal components of the magnetic field according to (|8jl. 

(A18) 

(A19) 

Both components of the field vary with radius as r u ~ 2 , as 
they should for the assumed form of P(r,9). Moreover, B r 
does not vary with 9, which i s approximately the case for 
self-similar force-free models l|Naravan et al.lT2007h . In the 
context of our problem, we see that the normal component of 
the field at the surface of the central compact object (r = 1) 
is independent of 9. Thus we have a uniform field emerging 
from the object — a split monopole. This is a nice property 
of the flux function (|A8[) . 

For reference, we provide the poloidal components of 
the magnetic field in cylindrical coordinates, 

q 

Br — B r sin 8 + Bg cos 9 = r"~ [1 + (1 — v) cos 6] tan — , 

Q 

B z — B r cos 9 — Bg sin 9 — r v ~ [cos 9 + v sin Stan — ], 

In the disc plane (9 = 7r/2), we have 

Br = r v ~ 2 = B r , (A20) 

B z = vr v ~ 2 = -Bg. (A21) 

The normal component of the field on the disc varies as a 
power-law in radius. 

A4 Lorentz Factor 

A4-1 General Formula for the Lorentz Factor 

For an axisymmetric configuration in steady state, the angu- 
lar velocity Q is constant on field lines. This means that the 
comoving frame of the fluid at any position rotates at the 
same fi as the angular velocity of the disc at the correspond- 
ing footpoint. In this rotating frame, there is no electric field 
(since we assume infinite conductivity). Therefore, in the in- 
ertial/lab frame, the electric fiel d E is perpendicular to the 
magnetic field B and is eq ual to l|Thorne et al.| [l986; Bcskin 
ll997l ; lNaravan et alj|2007f ) 



E = S1RB V 



(A22) 



where B p is the magnitude of the poloidal magnetic field. 

In the force- free approximation, the velocity of the 
fluid may be conven iently taken to be the drift velocity 
l|Naravan et alj|2007| l. 



Vd 



E x B 



B 2 



E 
B 



This gives the Lorentz factor 7 of the fluid, 



1 



B 



E 2 



B- 



Bp 
B 2 



+ 



B 



E l 



B 2 



(A23) 



(A24) 



Consider now the far asymptotic region of the jet where 
Q.R ^> 1 and the velocity is ultrarelativistic, s 3 fs 1, 7 > 1. 
In this limit we have E S> B p due to (|A22|) and E « B due 
to (|A23|) . These two relations, combined with the defini- 
tion of B in terms of B v and B v , allow us to obtain the 
relative magnitudes of the electromagnetic field components, 



B^\B V \^E^B P for QR > 1, v 2 



(A25) 



At small distances from the compact object the first term 
on the r.h.s. of (|A24[I dominates, and we write its inverse as 



B 2 



1 + 



U 2 

P 



7i 2 , 



where from (|A22|I and (|A25|) we have 

n2i 1/2 



71 



[1 + {my 



(A26) 



(A27) 



In this near-zone, we have 7 ~ 71 ~ £IR. We refer to 
this as the first acceleration regime. Formula (|A27|I . as we 
will see below, gives a good approximation to the Lorentz 
factor at moderate distances fro m the compact object and 
has been noted by many authors llBlandford fc Znaieklll977l; 
Beskin et aDl 19981 : iBeskin fc Nokhrinall2006l : iNaravan et all 
20071 ). 

For certain field geometries it is the second term on the 
r.h.s. of (|A24[l that determines the bulk Lorentz factor. The 
inverse of the second term can be written as 



B 



B 2 -E 2 



2 

= 72, 



(A28) 



where 72 is determined by the ratio of the local poloidal 
radius of curvature of the field line R c to t he cylindrical 
radius from the jet axis R l|Beskin et alj|2004 and CAT)) . 



72 w C 



1/2 



(A29) 



Here C is a numerical factor of order unity which has to 
be determined by a more careful study of the solution. The 
value of this factor depends on the spatial distribution of 
field line rotation frequency Q and as we show analytically 
in Appendix lA7l in the jet region, where £1 = const, we have 
C — v3 w 1.73. We find that this single value of C does a 
very good job of explaining all our numerical results, giving 
an error of less than 10% in all cases we have considered 
(see Sj3j . Equation ()A29|) corresponds to the second acceler- 
ation regime. It illustrates that, in a relativistic electromag- 
netic flow, collimation and acceleration are suppressed: for 
the fluid to have a sufficiently large Lorentz factor, the field 
line to which it is attached has to be sufficiently straight, 
i.e. have a sufficiently large value of R c /R. 

Combining equations ()A24|) - <|A29|I we obtain the fol- 
lowing general formula for the Lorentz factor, 



1 1 

7f 7 2 



(A30) 



where 71 and 72 are given by (|A27[) and (|A29[) . Although 
we derived this expression in the asymptotic limit of large 
Lorentz factors, the formula is quite accurate all the way 
from the foot-point of field lines to large distances (see J3j . 

A4-2 Lorentz Factor of the Jet and the Wind 

In order to numerically evaluate formula (|A30[) for our jet- 
wind solution, we adopt the field line shape QAllfl . In the 
jet all field lines rotate at the same frequency, 

n jct = Q.q = const, (A31) 

whereas in the wind field lines rotate differentially, according 
to the rotation profile in the disc, 



fT md (r,0) = n[r ip (r,9)] « Q o (r v 2 /2) 



"ft 2 l r >\~P/ l/ 



(A32) 
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We have used (|A17|I to substitute for ri p (r,8) and equa- 
tion to substitute for n(Rf p ). Combining expres- 
sions (IA31[) and (|A32|) we have 



n(r,0) 



no, 

^^o■ I '<9 2 /2)-' 3/ ^ 



(jet), 



" 3 "~" J \ (A33) 
> Oj (wind). 



The first acceleration regime is straightforward. For 
both the jet and the wind we have 



71 = Q(r,0)rO. 



(A34) 



For the second acceleration regime, we use the asymptotic 
field line shape <)A13[) to compute the curvature radius of 
the field line, 



Rc 



[1 + (dR/dz 



,21 3/2 



\d 2 R/dz 2 \ \d?R/dz 2 Y 
which together with the field line shape (|A13f) gives 



Rc 
rip 



\ l+u/2 



r {p J 



(1- t//2)i/sin(0fe/2) 



(A35) 



(A36) 



Therefore, according to (|A13I) and (|A29[I . for the jet field 
lines we have 



72 



2C 



y/(2-u)u' 



(A37) 



Note that the only explicit dependence of the Lorentz factor 
on position in the second acceleration regime is the spherical 
polar angle. Even though this expression rather accurately 
describes the distribution of the Lorentz factor in the jet 
with one single value of C = v3 (?JO}> we find that in the 
wind the values of C are different for different field lines 
because in the wind the field line shape deviates from the 
analytical expectation (|A13|l . 



A4-3 Lorentz Factor Scaling Along Field Lines 

Consider a field line starting at a foot-point r = n p , 8 = 
6{p, and rotating at a frequency S7f p . For a jet field line we 
have Tf p = 1 whereas for a wind field line we have 8f p = 
7r/2. Let us evaluate the Lorentz factor scalings (|A27|I and 
(TA291) along a given field line as a function of distance from 
the compact object. Using various results derived earlier we 
obtain 



71 « fifp n p I — ) 
\ripj 



l-v/2 



x 2 sin 



9ft, 



,/2 



(A38) 
(A39) 



Since C = v3 and flf p is substantially less than unity, 
it is easy to see that at small radii we have 71 < 72. By 
equation (|A30[) . 7 is determined by the smaller of 71 and 
72. Therefore, near the base of each field line the fluid is 
always in the first acceleration regime. As the field line goes 
to larger distances, 71 and 72 increase at different rates. If 
v ^ I, 71 always remains smaller than 72, and the first accel- 
eration regime operates throughout the field line. However, 
for v < 1, 72 grows more slowly than 71. At a certain dis- 
tance along the field line, 72 then becomes smaller than 71 
and the flow makes a transition to the second acceleration 



regime (jNaravan et alj l2007t ) . The transition between the 
two regimes happens at a distance rti given by 



rtr 
rip 



C 



n fp rt p 2sin 2 (0 fp /2)y(2 



1/(1-") 



(A40) 



Equation (|A40|) shows that the transition occurs sooner 
for field lines with foot-points at lower latitudes (9f p closer 
to 7r/2) on the compact object compared to field lines that 
emerge near the pole. This leads to the development of a 
fast jet core in the polar region, where the Lorentz factor 
continues to be determined by the first acceleration regime 
for a long distance, surrounded by a slower sheath, where 
the acceleration switches to the second less-efficient regime 
at a shorter distance (see Sj3|. 



A5 Toroidal Magnetic Field and Enclosed 
Poloidal Current 

To obtain the toroidal magnetic field B^, we use the fact 
that in ideal MHD the enclosed poloidal current, 

/ = RB v /2, (A4I) 

is constant along each fiel d line (|Thorne et al.|[l986l ; iBeskinl 
ll997l ; lNaravan et al.ll2007T l. A positive value of I indicates a 
current in the positive z direction. 

Asymptotically far from the compact object, where the 
fluid motion is ultrarelativistic (7 > I), we have according 
to (TA25|) . 



B - 



-QRB V 



(A42) 



where the negative sign is because the field lines are swept 
back in the opposite direction to rotation, i.e. in the negative 
(^-direction for a positive Q. Since asymptotically, 8 <C 1, 
according to (|AI8|) - (|AI9|) we have B p w B r = r"~ 2 and 
therefore 



72 



-fi(r"(9 2 /2), 



(A43) 



where we have approximated R ss r8. According to (|A8)> . 
the expression in parentheses is approximately equal to the 
flux function P. Therefore we have 



-np. 



(A44) 



This relation between two quantities that are each con- 
served along field lines must hold throughout the solution 
(although we derived it only asymptotically). Therefore, ac- 
cording to (|A41|) and (|A44[) . we obtain 



2np 



2JV(1- 



2 / 
(A45) 



A6 Magnetic Pressure 

Because of the assumption of perfect conductivity, the elec- 
tric field in the comoving frame of the fluid vanishes and 
there is only a comoving magnetic field b. Given the electric 
and magnetic fields in the lab frame, E and B, the comoving 
magnetic field strength is 



b = ^B 2 - E 2 , 



(A46) 
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and the comoving magnetic pressure is 



From (|A46|) . we can write 



6 2 = + s! 



£ 2 



R 2 B 

P ' 2 ' 
72 



(A47) 



(A48) 



where the approximate equality is due to (|A28|I , Further, 
due to (|A25|I and (|A27|) . we can rewrite (|A47|) as 



(A49) 



where B 2 / (8w) is the pressure that the jet would have if the 
compact object was not spinning (when B v = E = 0). 

If v ^ 1, we saw earlier that we always have 71 < 72. 
Therefore, using the expression for the radial magnetic 
field (|A18|) . we obtain 



v > 1 



Pmag OC T 



2(v-2) 



(A50) 



The magnetic pressure follows a single power-law. However, 
for v < 1, the term involving the Lorentz factors in (|A48|) 
cannot be neglected. Using the expressions for 71 (|A38| . 
72 (|A39[) and the transition radius r tr (IA40|) . we obtain 



2(l-f) 



(A51) 



and the magnetic pressure behaves as a broken power-law, 
^ < 1 : Pmae oc ^ ' „ (A52) 



r > 



Note that the power-law index at large distances is indepen- 
dent of the value of v. 



A7 Transverse Force Balance 

In this section we study the steady state transversal force 
balance in force-free jets and analytically determine the 
value of the prefactor C that scales the Lorentz factor in 
the second acceleration regime (eq. I A29|) . The force balance 
equation is 



j x B + pE = 0, 



(A53) 



where p and j are the lab-frame electric charge and current 
density. Using Ampere's law, 



(A54) 



we have 



iwj xB = CVxB)xB — (BV)B - - V(B ). (A55) 
Let us introduce a local orthonormal tetrad 

(t = B p /B p ,0,n = E/E), (A56) 
corresponding to local rotated cylindrical coordinates, 

(T,Rp,n), (A57) 



where r measures the distance along field lines, n measures 
distance perpendicular to field lines, and Rip measures the 
distance in the toroidal direction. We have then 

(BV)B = (Br-jL + B v -£- + B n ^-){B T f + B^ip + B n h). 
or Roip on 

Projecting the above relation along h, we obtain: 

n ■ [(BV)B] = n ■ B T ^-{B T f) + h ■ B v ^-{B v 0) 



(A58) 



where R is a unit vector directed along the cylindrical radius. 
Combining the above results, we have for the force balance 
in the n-direction: 



AtyR c 4ttR ( ' dn\8TT 



47T 



E = 0. (A59) 



We can write the divergence of the electric field in the local 
rotated cylindrical coordinates as 



dE L + }_dR E | dE„ 
dr R dn ™ dn ' 



-sr. (A60) 



where we have used the fact that E T — E v — 0. The first 
term in (]A60[) can be rewritten in terms of the curvature 
radius of the field line: 



dE T 
dr 



d{Eh) _ / dh 
dr \ dr 



E_ 
R~c 



(A61) 



Finally, combining equations (|A59|I - (|A61|I . we obtain 
a form of the force balance equation that does not explicitly 
contain any derivatives along field lines: 



Bp - E 1 
4ttR c 



R 2 

(R-h)^ 



E 2 



A-kR 



d_ 

dh 



B l - E' 
8^ 



0. (A62) 



This equation is the analog of equation (79) in 
Bcs kin fc Nokhrinal (120061 . they have a sign typo) . 

At large distances from the compact object we have 
(R ■ h) » -1 and d/dh » -d/dR (ea.\Mty. With this we 
can re-write (|A62|) in projection along R, 



B2 T?2 r>2 rriz 



4tvR c 



AtvR 



d_ 

dR 



Bp + B ¥ 



E 2 



871 



0. (A63) 



This form explicitly shows that the poloidal and toroidal 
relativistic hoop stresses, which are the first two terms, are 
balanced by the comoving pressure gradient, the last term. 

Assuming that the field line rotation frequency varies 
as a power-law in the cylindrical radius, Q oc i? A , we have 
B v w —QRBp oc R 1+x (eq.EHJ, where we took Bp to be 
nearly independent of R (eg. IA18I - IA19|) . With the help 
of (|A28|) we have, 



B, 



oc R 



A+2\ 



(A64) 



where we used eq. (|A37[) for the Lorentz factor in the second 
acceleration regime, which gives 72 oc 1/R. We note that in 
the wind the poloidal magnetic field is not uniform, and 72 
follows a different power-law. 



© 2008 RAS, MNRAS 000. [111251 



Simulations of Ultrarelativistic Magnetodynamic GRB Jets 23 



In the second acceleration regime we have 71 3> 72 and 
thus B 2 -E 2 > B 2 (Appendix lA4.1| l. Therefore in this limit 
the terms involving B p in (|A63[) can be dropped, 



B; 



,-E 2 d_ 
R + dR 



E* 



F 2 

— . (A65) 



By combining equations (|A64fl and (|A65|) with (IA25|) 
and (fA28l . we can express the Lorentz factor in the second 
acceleration regime as 



72 



E 2 \ 1/2 (R^ 1 ' 2 
B&-&) * \ R 



where 



c = Vm, A : 



dlog» 
dlogR' 



(A66) 



(A67) 



Specializing to the jet region where f2 = const, or A = 0, we 
obtain C = \/H ~ 1.7. Numerically, however, we find a range 
of values for C rather than just a single value ( i]3.3fl . One 
of the reasons for this may be the fact that instead of be- 
ing exactly constant in the jet, fl(R) slightly decreases with 
increasing R near the jet-wind boundary due to numerical 
diffusion (see Fig. [6j . This causes A to vary from at the jet 
axis to ~ — 1 near the jet-wind boundary (Fig. [6}. Accord- 
ing to (|A67I) . this variation in A translates into a variation 
in C = \/2 — %/3 ~ 1.4 — 1.7 across the jet, consistent with 
the numerical simulations ( i]3.3l Fig. II Of) . 



A8 Power Output 

The power output in a magnetodynamic jet is equal to the 
Poynting flux, whose magnitude is given by 



S = 



\E x B\ 
4tt 



(A68) 



At asymptotically large distances, according to (|A25[l we 
can approximately write 

B ^E^ flRB p . (A69) 

Therefore, using (lA33]l . we have 

n(r,9) 2 R 2 B 2 {l(r,8) 2 r 2 »- 2 9 2 
4^ 4^ ' (A70) 

where we have approximated B p ~ B r ~ r v ~~ 2 with the help 
of (|A18I) . We can now evaluate the lateral dependence of S 
in the jet and the wind, 



2J2v-2n2 



n 2 r 



1- \nlr 2v - 2 e 2 {r v e 2 /2y 



2/3/v 



> 



(A71) 



Let us define the energy flux per unit solid angle as 



dP 

dui 



r 2 S{r, 0). 



(A72) 



For small angles the energy flux grows quadratically with the 
polar angle inside the jet, tx 8 2 , and then falls off rapidly in 
the wind, tx g 2-4 ' 3 /" _ The maximum power output occurs at 
the jet- wind boundary, — 0j, 



dP 

du) 



(A73) 



These results are confirmed in Fig. 1111 Keplerian discs have 



— —3/2, so the power output in the wind falls off ex- 
tremely rapidly (e.g. tx 6~ 6 for our fiducial model A which 
has v = 0.75). However, the contribution to the total energy 
output from the disc may still be significant. 
The total power output of the jet is 



P jeI (r) = 



dP 

2tt sin6» dO i 

doj 



1 

32' 



n 2 



r v 6< 



(A74) 



The region 6 < 6j = \p2r~ v l 2 defines the jet region, 
and the numerical evaluation was performed for the fidu- 
cial model A. For the disc wind, 



(r) 



dP 

2-Ksme — de i 
doj 



n 2 



r fQ2\ 2-2/3/" 



or, 



209/iz -1) 



1_ 

32' 



(A75) 



where the numerical evaluation, again, was performed for 
model A. For this model the disc wind provides the same 
energy output as the jet. 

Converting the results to physical units, the total power 
output of a jet in an astrophysical system with BH mass M, 
radial magnetic field strength near the BH B r , and angular 
rotation frequency fio, is 



P jct = -tf B, 
2 



2 2 
r c 



1.8 x 10° 



( M 



10 15 G J \ 3M F 



[?]■ 

(A76) 



where ro 



= GM/c 2 for a rapidly spinning BH, and 



f2 max = 0.25c/r 9 is the maximum angular rotation frequency 
that magnetic field lines threading a spinning compact ob- 
ject can have in a stationary state. 



APPENDIX B: FLUID SPEED IN FORCE-FREE 
EQUILIBRIUM 

In this section we would like to answer the following ques- 
tion: In the limit of infinite magnetization, i.e., a — > 00, are 
we guaranteed that a force-free solution will give a physi- 
cally meaningful solution for any fluid that is carried along 
with the field? Specifically, if we have a force-free solution 
which has drift speed Vd = E/B feq, IA23]l less than c every- 
where of interest, and if we add some fluid with negligible 
inertia [a — > 00) which is carried along with the field, will 
the flow speed Vf of the fluid be physical, i.e., will we have 
Vf < c everywhere? The answer to this question is not obvi- 
ous. A fluid possesses additional conserved quantities, e.g., 
the Bernoulli constant, which are not relevant for the elec- 
tromagnetic field. Could the additional constraints give in- 
consistent results for the fluid motion even though the fields 
behave physically? 

Consider a steady axisymmetric force-free solution in 
which the poloidal field is B p h p along a poloidal unit vector 
n p and the toroidal field is B^h^ along the toroidal unit 
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vector fitp. Assume that B p is positive, and write the total 
field strength as 



B = (Bp + B 



2\l/2 



aB v , a > 1. 



(Bl) 



Assume further that the rotation is in the positive sense, 
which means that the toroidal component of the field will 
be negative: 



B v = —(B — Bp) ' = 



20/2 _ (a 2 ~ 1) 



1/2 



-B = -{a 2 



a/2 



Br, 



(B2) 

Let us write the comoving magnetic field in terms of the 
lab-frame magnetic field, 

b — (B 2 — E 2 ) 1 ^ 2 = SB, 0<5^1. (B3) 

We will see below that S determines the drift velocity of 
magnetic field. The electric field is given by 



SIR 



SIRB V 



(B4) 
(B5) 



From this it follows that 

SIR = a(l - S 2 ) 1/2 c. 
We then obtain the drift speed, 

E x B SIRB P , „ n A , 

-— — {—B p n v + B v n p ) = Vdpnp+Vd^ntp, 

(B6) 

where the two velocity components are given by 
(a 2 -l)Va(i_jaji/a 



v d = c- 



B 1 



Vdp = 

(1-52)1/2 
Vdip = C. 



(B7) 
(B8) 



The Lorentz factor corresponding to the drift speed is then 

f f 



7 = 



(B9) 



(l-«2/£?)Va 5- 

We now see that 5 is simply the inverse of the drift Lorentz 
factor. 

In ^A4.1 1 we showed that there are two distinct regimes 
of acceleration in the jet. In the first regime, we have 7 = 
(I + Sl 2 R 2 /c 2 ) 1/2 . Making use of equations Q35) and (|B9|) . 
we obtain 



First regime : 5 = — . 



(BIO) 



In the second regime, the Lorentz factor depends on the 
shape of the field line, i.e., it is no longer a local quantity 
but depends on spatial derivatives of the field components. 
Nevertheless, we know that 7 in this regime is smaller than 
for the first regime. Therefore, we have 

Second regime : 5 > — . (Bff) 



Combining the two, we have the condition 
Either regime : 5 ^ —. 



(B12) 



The fluid velocity Vf is constrained to be equal to 
SIR hip plus an arbitrary velocity parallel to B. This is 
expressed by the following condition on the poloidal and 
toroidal components of the fluid velocity 

v Sv = SlR+^ Vfp = a(l-S 2 ) 1/2 c-(a 2 ~l) 1/2 v fp . (B13) 

Br, 



Note in particular that the two components of the drift ve- 
locity satisfy this condition. Let us now rewrite the fluid 
velocity in terms of the drift velocity as follows, 



Vfp = v dp + <lc = 

"(!_ £2)1/2 



(a 2 -l) 1/2 (l-<5 2 ) 



2U/2 



+ t 



-(a 2 -I 



,1/2, 



(BI4) 
(B15) 



where e is a small parameter which determines the deviation 
of the fluid velocity from the drift velocity. The fluid Lorentz 
factor is then 



7/ 



1 



1 



il-(v% + v 2 f ^)/c 2 ] 1/2 (S 2 



2)1/2 



(B16) 



Thus, if we want the fluid motion to be physical (vf < c), 
the following condition must be satisfied, 



e < 



1 

0:7 



(B17) 



We now use the Bernoulli equation to fix the 
value of the parameter e. According to the Bernoulli 
and angular momentum conservation equations, the fol- 
lowin g quantity must be con s erved along each field 
line (Bckcnste in fc Oronl 1 19781: iMestel fc Shibatal 1 19941 : 
IContopouloslll995l : IContopoulos et al1ll999h : 



M = 7/ 



1 - 



SlRvf v 



= constant. 



(B18) 



Near the base of the jet the flows that we consider in this 
paper are sub-Alfvenic, so we expect the fluid motion to be 
not very relativistic. Therefore, we expect fi to be ~ 1. For 
instance, if the fluid moves at the drift velocity at the base 
of the jet, then we have 



M = 



7/ 



SlRvf^ 



So = — , 

70 



(B19) 



where the subscript refers to values at the base of the jet. 
For most problems of interest to us, 70 ~ 1, so we have 
A*~ I- 

Now let us apply the jj, constraint at a general point in 
the force-free jet. This gives the condition 



M = 



5 2 + a(l-a 2 ) 1/2 (a 2 -1) 1/2 

(52_ a 2 £ 2)l/2 



(B20) 



which is a quadratic equation for e. The quadratic is easily 
solved and gives the somewhat messy result, 

= s(» 2 - s 2 ) 

6 Q^(Q 2 + /i 2 - 1 - a 2 5 2)l/2 + q5 ( q 2 _ 1)1/2(1 _ £2)1/2 ' 

(B21) 

So far, we have not made any approximations; all the 
results are exact. Now, for simplicity, let us assume that 
we are in an asymptotic region of the jet where a 2> 1 
and 7 3> 1. In this limit, the first term in the denominator 
of (|B21I) dominates over the second. We thus have 



J±_ 

q 2 7 



(B22) 



which clearly satisfies the condition (|B17|1 that implies vj < 
c. Thus if Vd < c in the force-free solution, then the total 
fluid velocity satisfies Vf < c even in the limit of a — > 00. 
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Substituting this solution into ()B 16|) . we find that the devia- 
tion of the fluid Lorentz factor from the drift Lorentz factor 
is 

A7 = 7/ _ 7 w _L_ f <<c L (B23) 

The inequality on the right follows from the fact that a ^ 
1/5 = 7 > 1 (see eq. IB12f) and fi ~ 1. 

Equation f|B23f) shows that any fluid which is carried 
along with an infinitely magnetized force-free flow has only 
a slightly modified Lorentz factor relative to the drift speed. 
In other words, the perturbation to the plasma motion is 
negligibly small and we can, for all practical purposes, as- 
sume that the fluid has the same speed as the drift speed of 
the electromagnetic fields. 
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